ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ComponentGGNuclNuclXsc.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ComponentGGNuclNuclXsc.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 // 24.11.08 V. Grichine - first implementation
27 //
28 // 04.09.18 V. Ivantchenko Major revision of interfaces and implementation
29 // 27.05.19 V. Ivantchenko Removed obsolete methods and members
30 
32 
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4NucleiProperties.hh"
36 #include "G4ParticleDefinition.hh"
37 #include "G4HadronNucleonXsc.hh"
39 #include "G4NuclearRadii.hh"
40 
41 static const G4double inve = 1./CLHEP::eplus;
42 
44  : G4VComponentCrossSection("Glauber-Gribov Nucl-nucl"),
45  fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0), fProductionXsc(0.0),
46  fDiffractionXsc(0.0), fEnergy(0.0), fParticle(nullptr), fZ(0), fA(0)
47 {
50  fHNXsc = new G4HadronNucleonXsc();
52 }
53 
55 {
56  delete fHNXsc;
57 }
58 
60 
62  const G4ParticleDefinition* aParticle, G4double kinEnergy,
63  G4int Z, G4double A)
64 {
65  ComputeCrossSections(aParticle, kinEnergy, Z, G4lrint(A));
66  return fTotalXsc;
67 }
68 
70 
72  const G4ParticleDefinition* aParticle, G4double kinEnergy,
73  G4int Z, G4int A)
74 {
75  ComputeCrossSections(aParticle, kinEnergy, Z, A);
76  return fTotalXsc;
77 }
78 
80 
82  const G4ParticleDefinition* aParticle, G4double kinEnergy,
83  G4int Z, G4double A)
84 {
85  ComputeCrossSections(aParticle, kinEnergy, Z, G4lrint(A));
86  return fInelasticXsc;
87 }
88 
90 
92  const G4ParticleDefinition* aParticle, G4double kinEnergy,
93  G4int Z, G4int A)
94 {
95  ComputeCrossSections(aParticle, kinEnergy, Z, A);
96  return fInelasticXsc;
97 }
98 
100 
102  const G4ParticleDefinition* aParticle, G4double kinEnergy,
103  G4int Z, G4double A)
104 {
105  ComputeCrossSections(aParticle, kinEnergy, Z, G4lrint(A));
106  return fElasticXsc;
107 }
108 
110 
112  const G4ParticleDefinition* aParticle, G4double kinEnergy,
113  G4int Z, G4int A)
114 {
115  ComputeCrossSections(aParticle, kinEnergy, Z, A);
116  return fElasticXsc;
117 }
118 
120 
122  const G4ParticleDefinition* aParticle, G4double kinEnergy,
123  G4int Z, G4int A)
124 {
125  ComputeCrossSections(aParticle, kinEnergy, Z, A);
126  return (fInelasticXsc > fProductionXsc)
128 }
129 
131 
133 {}
134 
136 
138 {
139  G4cout << "G4ComponentGGNuclNuclXsc: uses Glauber-Gribov formula" << G4endl;
140 }
141 
143 
144 void G4ComponentGGNuclNuclXsc::Description(std::ostream& outFile) const
145 {
146  outFile << "G4ComponentGGNuclNuclXsc calculates total, inelastic and\n"
147  << "elastic cross sections for nucleus-nucleus collisions using\n"
148  << "the Glauber model with Gribov corrections. It is valid for\n"
149  << "all incident energies above 100 keV./n"
150  << "For the hydrogen target G4HadronNucleonXsc class is used.\n";
151 }
152 
154 //
155 // Calculates total and inelastic Xsc, derives elastic as total - inelastic
156 // accordong to Glauber model with Gribov correction calculated in the dipole
157 // approximation on light cone. Gaussian density of point-like nucleons helps
158 // to calculate rest integrals of the model. [1] B.Z. Kopeliovich,
159 // nucl-th/0306044 + simplification above
160 
162  const G4ParticleDefinition* aParticle, G4double kinEnergy,
163  G4int Z, G4int A)
164 {
165  // check cache
166  if(aParticle == fParticle && fZ == Z && fA == A && kinEnergy == fEnergy)
167  { return; }
168  fParticle = aParticle;
169  fZ = Z;
170  fA = A;
171  fEnergy = kinEnergy;
172 
173  G4int pZ = G4lrint(aParticle->GetPDGCharge()*inve);
174  G4int pA = aParticle->GetBaryonNumber();
175 
176  // hydrogen
177  if(1 == Z && 1 == A) {
178  G4double e = kinEnergy*CLHEP::proton_mass_c2/aParticle->GetPDGMass();
185  return;
186  }
187  static const G4double cofInelastic = 2.4;
188  static const G4double cofTotal = 2.0;
189 
190  G4double pTkin = kinEnergy/(G4double)pA;
191 
192  G4int pN = pA - pZ;
193  G4int tN = A - Z;
194 
195  G4double tR = G4NuclearRadii::Radius(Z, A);
196  G4double pR = G4NuclearRadii::Radius(pZ, pA);
197 
198  G4double cB = ComputeCoulombBarier(aParticle, kinEnergy, Z, A, pR, tR);
199 
200  if ( cB > 0. )
201  {
202  G4double sigma = (pZ*Z+pN*tN)*fHNXsc->HadronNucleonXscNS(theProton, theProton, pTkin);
204 
205  sigma += (pZ*tN+pN*Z)*fHNXsc->HadronNucleonXscNS(theNeutron, theProton, pTkin);
207 
208  // G4cout<<"ppInXsc = "<<ppInXsc/millibarn<<"; npInXsc = "<<npInXsc/millibarn<<G4endl;
209  // G4cout<<"npTotXsc = "<<fHNXsc->GetTotalHadronNucleonXsc()/millibarn<<"; npElXsc = "
210  // <<fHNXsc->GetElasticHadronNucleonXsc()/millibarn<<G4endl;
211 
212  G4double nucleusSquare = cofTotal*CLHEP::pi*( pR*pR + tR*tR ); // basically 2piRR
213 
214  G4double ratio= sigma/nucleusSquare;
215  fTotalXsc = nucleusSquare*G4Log( 1. + ratio )*cB;
216  fInelasticXsc = nucleusSquare*G4Log( 1. + cofInelastic*ratio )*cB/cofInelastic;
218 
219  G4double difratio = ratio/(1.+ratio);
220  fDiffractionXsc = 0.5*nucleusSquare*( difratio - G4Log( 1. + difratio ) );
221 
222  G4double xratio= ((pZ*Z+pN*tN)*ppInXsc + (pZ*tN+pN*Z)*npInXsc)/nucleusSquare;
223  fProductionXsc = nucleusSquare*G4Log( 1. + cofInelastic*xratio)*cB/cofInelastic;
225  }
226  else
227  {
228  fInelasticXsc = 0.;
229  fTotalXsc = 0.;
230  fElasticXsc = 0.;
231  fProductionXsc = 0.;
232  fDiffractionXsc= 0.;
233  }
234 }
235 
237 
239  const G4ParticleDefinition* aParticle,
240  G4double pTkin, G4int Z, G4int A,
241  G4double pR, G4double tR)
242 {
243  G4int pZ = aParticle->GetPDGCharge()*inve;
244  G4double pM = aParticle->GetPDGMass();
246  G4double pElab = pTkin + pM;
247  G4double totEcm = std::sqrt(pM*pM + tM*tM + 2.*pElab*tM);
248  G4double totTcm = totEcm - pM -tM;
249 
250  static const G4double qfact = CLHEP::fine_structure_const*CLHEP::hbarc;
251  G4double bC = qfact*pZ*Z*0.5/(pR + tR);
252 
253  G4double ratio = (totTcm <= bC ) ? 0. : 1. - bC/totTcm;
254  // G4cout<<"G4ComponentGGNuclNuclXsc::ComputeCoulombBarier= "<<ratio
255  // <<"; pTkin(GeV)= " <<pTkin/GeV<<";
256  // " pPlab = "<<pPlab/GeV<<"; bC = "<<bC/GeV<<"; pTcm = "
257  // <<pTcm/GeV<<G4endl;
258  return ratio;
259 }
260 
262 //
263 // Return single-diffraction/inelastic cross-section ratio
264 
266  const G4DynamicParticle* aParticle, G4double tA, G4double tZ)
267 {
268  ComputeCrossSections(aParticle->GetDefinition(),
269  aParticle->GetKineticEnergy(),
270  G4lrint(tZ), G4lrint(tA));
271 
272  return (fInelasticXsc > 0.0) ? fDiffractionXsc/fInelasticXsc : 0.0;
273 }
274 
276 //
277 // Return quasi-elastic/inelastic cross-section ratio
278 
280  const G4DynamicParticle* aParticle, G4double tA, G4double tZ)
281 {
282  ComputeCrossSections(aParticle->GetDefinition(),
283  aParticle->GetKineticEnergy(),
284  G4lrint(tZ), G4lrint(tA));
285 
286  return (fInelasticXsc > 0.0) ? 1.0 - fProductionXsc/fInelasticXsc : 0.0;
287 }
288