| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // 24.11.08 V. Grichine - first implementation
27 //
28 // 04.09.18 V. Ivantchenko Major revision of i
29 // 27.05.19 V. Ivantchenko Removed obsolete me
30
31 #include "G4ComponentGGNuclNuclXsc.hh"
32
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4NucleiProperties.hh"
36 #include "G4ParticleDefinition.hh"
37 #include "G4HadronNucleonXsc.hh"
38 #include "G4ComponentGGHadronNucleusXsc.hh"
39 #include "G4NuclearRadii.hh"
40 #include "G4Pow.hh"
41
42 namespace
43 {
44 const G4double inve = 1./CLHEP::eplus;
45 }
46
47 G4ComponentGGNuclNuclXsc::G4ComponentGGNuclNuc
48 : G4VComponentCrossSection("Glauber-Gribov Nu
49 {
50 theProton = G4Proton::Proton();
51 theNeutron = G4Neutron::Neutron();
52 theLambda = G4Lambda::Lambda();
53 fHNXsc = new G4HadronNucleonXsc();
54 fHadrNucl = new G4ComponentGGHadronNucleusXs
55 }
56
57 G4ComponentGGNuclNuclXsc::~G4ComponentGGNuclNu
58 {
59 delete fHNXsc;
60 }
61
62 //////////////////////////////////////////////
63
64 G4double G4ComponentGGNuclNuclXsc::GetTotalEle
65 const G4ParticleDefinition* aParticle
66 G4int Z, G4double A)
67 {
68 ComputeCrossSections(aParticle, kinEnergy, Z
69 return fTotalXsc;
70 }
71
72 //////////////////////////////////////////////
73
74 G4double G4ComponentGGNuclNuclXsc::GetTotalIso
75 const G4ParticleDefinition* aParticle
76 G4int Z, G4int A)
77 {
78 ComputeCrossSections(aParticle, kinEnergy, Z
79 return fTotalXsc;
80 }
81
82 //////////////////////////////////////////////
83
84 G4double G4ComponentGGNuclNuclXsc::GetInelasti
85 const G4ParticleDefinition* aParticle
86 G4int Z, G4double A)
87 {
88 ComputeCrossSections(aParticle, kinEnergy, Z
89 return fInelasticXsc;
90 }
91
92 //////////////////////////////////////////////
93
94 G4double G4ComponentGGNuclNuclXsc::GetInelasti
95 const G4ParticleDefinition* aParticle
96 G4int Z, G4int A)
97 {
98 ComputeCrossSections(aParticle, kinEnergy, Z
99 return fInelasticXsc;
100 }
101
102 //////////////////////////////////////////////
103
104 G4double G4ComponentGGNuclNuclXsc::GetElasticE
105 const G4ParticleDefinition* aParticle
106 G4int Z, G4double A)
107 {
108 ComputeCrossSections(aParticle, kinEnergy, Z
109 return fElasticXsc;
110 }
111
112 //////////////////////////////////////////////
113
114 G4double G4ComponentGGNuclNuclXsc::GetElasticI
115 const G4ParticleDefinition* aParticle
116 G4int Z, G4int A)
117 {
118 ComputeCrossSections(aParticle, kinEnergy, Z
119 return fElasticXsc;
120 }
121
122 //////////////////////////////////////////////
123
124 G4double G4ComponentGGNuclNuclXsc::ComputeQuas
125 const G4ParticleDefinition* aParticle
126 G4int Z, G4int A)
127 {
128 ComputeCrossSections(aParticle, kinEnergy, Z
129 return (fInelasticXsc > fProductionXsc)
130 ? (fInelasticXsc - fProductionXsc)/fInelas
131 }
132
133 //////////////////////////////////////////////
134
135 void G4ComponentGGNuclNuclXsc::BuildPhysicsTab
136 {}
137
138 //////////////////////////////////////////////
139
140 void G4ComponentGGNuclNuclXsc::DumpPhysicsTabl
141 {
142 G4cout << "G4ComponentGGNuclNuclXsc: uses Gl
143 }
144
145 //////////////////////////////////////////////
146
147 void G4ComponentGGNuclNuclXsc::Description(std
148 {
149 outFile << "G4ComponentGGNuclNuclXsc calcula
150 << "elastic cross sections for nucle
151 << "the Glauber model with Gribov co
152 << "all incident energies above 100
153 << "For the hydrogen target G4HadronNucleo
154 }
155
156 //////////////////////////////////////////////
157 //
158 // Calculates total and inelastic Xsc, derives
159 // accordong to Glauber model with Gribov corr
160 // approximation on light cone. Gaussian densi
161 // to calculate rest integrals of the model. [
162 // nucl-th/0306044 + simplification above
163
164 void G4ComponentGGNuclNuclXsc::ComputeCrossSec
165 const G4ParticleDefinition* aParticle, G4
166 G4int Z, G4int A)
167 {
168 // check cache
169 if(aParticle == fParticle && fZ == Z && fA =
170 { return; }
171 fParticle = aParticle;
172 fZ = Z;
173 fA = A;
174 fEnergy = kinEnergy;
175 G4Pow* pG4Pow = G4Pow::GetInstance();
176
177 G4int pZ = G4lrint(aParticle->GetPDGCharge()
178 G4int pA = aParticle->GetBaryonNumber();
179 G4int pL = aParticle->GetNumberOfLambdasInHy
180 G4bool pHN = aParticle->IsHypernucleus();
181 G4double cHN(0.88);
182
183 // hydrogen
184 if(1 == Z && 1 == A) {
185 G4double e = kinEnergy*CLHEP::proton_mass_
186 fHadrNucl->ComputeCrossSections( theProton
187 fTotalXsc = fHadrNucl->GetTotalGlauberGrib
188 fElasticXsc = fHadrNucl->GetElasticGlauber
189 fInelasticXsc = fHadrNucl->GetInelasticGla
190 fProductionXsc = fHadrNucl->GetProductionG
191 fDiffractionXsc = fHadrNucl->GetDiffractio
192 return;
193 }
194 static const G4double cofInelastic = 2.4;
195 static const G4double cofTotal = 2.0;
196
197 G4double pTkin = kinEnergy/(G4double)pA;
198
199 G4int pN = pA - pZ;
200 G4int tN = A - Z;
201
202 G4double tR = G4NuclearRadii::Radius(Z, A);
203 G4double pR = G4NuclearRadii::Radius(pZ, pA)
204
205 if(pHN)
206 pR *= std::sqrt( pG4Pow->Z23( pA - pL ) +
207
208 G4double cB = ComputeCoulombBarier(aParticle
209
210 if ( cB > 0. )
211 {
212 G4double sigma = (pZ*Z+pN*tN)*fHNXsc->Hadr
213 if(pHN) sigma += pL*A*fHNXsc->HadronNucleo
214 G4double ppInXsc = fHNXsc->GetInelasticHad
215
216 sigma += (pZ*tN+pN*Z)*fHNXsc->HadronNucleo
217 G4double npInXsc = fHNXsc->GetInelasticHad
218
219 G4double nucleusSquare = cofTotal*CLHEP::p
220
221 G4double ratio= sigma/nucleusSquare;
222 fTotalXsc = nucleusSquare*G4Log( 1. +
223 fInelasticXsc = nucleusSquare*G4Log( 1. +
224 fElasticXsc = std::max(fTotalXsc - fInel
225
226 G4double difratio = ratio/(1.+ratio);
227 fDiffractionXsc = 0.5*nucleusSquare*( difr
228
229 G4double xratio= ((pZ*Z+pN*tN)*ppInXsc + (
230 fProductionXsc = nucleusSquare*G4Log( 1. +
231 fProductionXsc = std::min(fProductionXsc,
232 }
233 else
234 {
235 fInelasticXsc = fTotalXsc = fElasticXsc =
236 }
237 }
238
239 //////////////////////////////////////////////
240
241 G4double G4ComponentGGNuclNuclXsc::ComputeCoul
242 const G4ParticleDefinitio
243 G4double pTkin, G4int Z, G4int A,
244 G4double pR, G4double tR)
245 {
246 G4int pZ = aParticle->GetPDGCharge()*inve;
247 G4double pM = aParticle->GetPDGMass();
248 G4double tM = G4NucleiProperties::GetNuclear
249 G4double pElab = pTkin + pM;
250 G4double totEcm = std::sqrt(pM*pM + tM*tM +
251 G4double totTcm = totEcm - pM -tM;
252
253 // 0.5 defines shape of Cross section correc
254 // at cB = totTcm it become zero
255 static const G4double qfact = 0.5*CLHEP::elm
256 G4double bC = qfact*pZ*Z/(pR + tR);
257
258 G4double ratio = (totTcm <= bC) ? 0. : 1. -
259
260 #ifdef G4VERBOSE
261 if (GetVerboseLevel() > 1) {
262 G4cout << "G4ComponentGGNuclNuclXsc::Compu
263 << "; pTkin(GeV)=" << pTkin/CLHEP::MeV
264 << " totTcm= " << totTcm/CLHEP::MeV<< ";
265 << G4endl;
266 }
267 #endif
268 return ratio;
269 }
270
271 //////////////////////////////////////////////
272 //
273 // Return single-diffraction/inelastic cross-s
274
275 G4double G4ComponentGGNuclNuclXsc::GetRatioSD(
276 const G4DynamicParticle* aParticle, G
277 {
278 ComputeCrossSections(aParticle->GetDefinitio
279 aParticle->GetKineticEn
280 G4lrint(tZ), G4lrint(tA));
281
282 return (fInelasticXsc > 0.0) ? fDiffractionX
283 }
284
285 //////////////////////////////////////////////
286 //
287 // Return quasi-elastic/inelastic cross-sectio
288
289 G4double G4ComponentGGNuclNuclXsc::GetRatioQE(
290 const G4DynamicParticle* aParticle, G
291 {
292 ComputeCrossSections(aParticle->GetDefinitio
293 aParticle->GetKineticEn
294 G4lrint(tZ), G4lrint(tA));
295
296 return (fInelasticXsc > 0.0) ? 1.0 - fProduc
297 }
298
299 //////////////////////////////////////////////
300