ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NuclearRadii.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4NuclearRadii.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 // Geant4 class G4NuclearRadii
28 //
29 // Author V.Ivanchenko 27.05.2019
30 //
31 //
32 
33 #include "G4NuclearRadii.hh"
34 #include "G4Pow.hh"
35 #include "G4PhysicalConstants.hh"
36 #include "G4ParticleDefinition.hh"
37 #include "G4NucleiProperties.hh"
38 
39 G4Pow* G4NuclearRadii::fG4pow = G4Pow::GetInstance();
40 const G4double fAlpha = 0.5*CLHEP::fine_structure_const*CLHEP::hbarc;
41 const G4double fInvep = 1.0/CLHEP::eplus;
42 
43 G4double G4NuclearRadii::ExplicitRadius(G4int Z, G4int A)
44 {
45  G4double R = 0.0;
46  // Special rms radii for light nucleii
47  if(Z <= 4) {
48  if(A == 1) { R = 0.895*CLHEP::fermi; }// p
49  else if(A == 2) { R = 2.13*CLHEP::fermi; }// d
50  else if(Z == 1 && A == 3) { R = 1.80*CLHEP::fermi; }// t
51  else if(Z == 2 && A == 3) { R = 1.96*CLHEP::fermi; }// He3
52  else if(Z == 2 && A == 4) { R = 1.68*CLHEP::fermi; }// He4
53  else if(Z == 3) { R = 2.40*CLHEP::fermi; }// Li7
54  else if(Z == 4) { R = 2.51*CLHEP::fermi; }// Be9
55  }
56  return R;
57 }
58 
59 G4double G4NuclearRadii::Radius(G4int Z, G4int A)
60 {
61  G4double R = ExplicitRadius(Z, A);
62  if(0.0 == R) {
63  if (A <= 50) {
64  G4double y = 1.1;
65  if( A <= 15) { y = 1.26; }
66  else if( A <= 20) { y = 1.19; }
67  else if( A <= 30) { y = 1.12; }
68  G4double x = fG4pow->Z13(A);
69  R = y*(x - 1./x);
70  } else {
71  R = fG4pow->powZ(A, 0.27);
72  }
73  R *= CLHEP::fermi;
74  }
75  return R;
76 }
77 
78 G4double G4NuclearRadii::RadiusRMS(G4int Z, G4int A)
79 {
80  G4double R = ExplicitRadius(Z, A);
81  if(0.0 == R) {
82  R = 1.24*fG4pow->powZ(A, 0.28)*CLHEP::fermi;
83  }
84  return R;
85 }
86 
87 G4double G4NuclearRadii::RadiusNNGG(G4int Z, G4int A)
88 {
89  G4double R = ExplicitRadius(Z, A);
90  if(0.0 == R) {
91  if(A > 20) {
92  R = 1.08*fG4pow->Z13(A)*(0.85 + 0.15*G4Exp(-(G4double)(A - 21)/40.));
93  } else {
94  R = 1.08*fG4pow->Z13(A)*(1.0 + 0.3*G4Exp(-(G4double)(A - 21)/10.));
95  }
96  R *= CLHEP::fermi;
97  }
98  return R;
99 }
100 
101 G4double G4NuclearRadii::RadiusECS(G4int Z, G4int A)
102 {
103  G4double R=0.;
104  const G4double c[3]={0.77329745, 1.38206072, 30.28295235};
105  const G4double c1=c[0];
106  const G4double c2=c[1];
107  const G4double c3=c[2];
108 
109  // Special rms radii for light nuclei
110  if (A <= 30) {
111  G4double vn = 0.5*A + fG4pow->powN(0.028*A,2) - fG4pow->powN(0.011*A,3);
112  G4double dev = vn - (A-Z);
113  R = c1*fG4pow->Z13(A) + c2/fG4pow->Z13(A) + c3*dev*dev/(A*A);
114  } else if (A<=50){
115  G4double y = 1.1;
116  G4double x = fG4pow->Z13(A);
117  R = y*(x - 1./x);
118  }
119  return R*CLHEP::fermi;
120 }
121 
122 G4double G4NuclearRadii::RadiusHNGG(G4int A)
123 {
124  G4double R = CLHEP::fermi;
125  if(A > 20) {
126  R *= 1.08*fG4pow->Z13(A)*(0.8 + 0.2*G4Exp(-(G4double)(A - 20)/20.));
127  } else {
128  R *= 1.08*fG4pow->Z13(A)*(1.0 + 0.1*G4Exp(-(G4double)(A - 20)/20.));
129  }
130  return R;
131 }
132 
133 G4double G4NuclearRadii::RadiusKNGG(G4int A)
134 {
135  return 1.3*CLHEP::fermi*fG4pow->Z13(A);
136 }
137 
138 G4double G4NuclearRadii::RadiusND(G4int A)
139 {
140  G4double R = CLHEP::fermi;
141  if(1 == A) { return R*0.895; }
142  G4double x = R*fG4pow->Z13(A);
143  if(A <= 3.) { x *= 0.8; }
144  else { x *= 1.7; }
145  return R;
146 }
147 
148 G4double G4NuclearRadii::RadiusCB(G4int Z, G4int A)
149 {
150  G4double R = ExplicitRadius(Z, A);
151  if(0.0 == R) {
152  G4int z = std::min(Z, 92);
153  R = r0[z]*fG4pow->Z13(A)*CLHEP::fermi;
154  }
155  return R;
156 }
157 
158 G4double G4NuclearRadii::ParticleRadius(const G4ParticleDefinition* p)
159 {
160  G4double R = CLHEP::fermi;
161  G4int pdg = std::abs(p->GetPDGEncoding());
162  if(pdg == 2112 || pdg == 2212) { R *= 0.895; }
163  else if(pdg == 211) { R *= 0.663; }
164  else if(pdg == 321) { R *= 0.340; }
165  else { R *= 0.5; }
166  return R;
167 }
168 
169 G4double G4NuclearRadii::CoulombFactor(
170  const G4ParticleDefinition* theParticle,
171  const G4ParticleDefinition* nucleon,
172  G4double ekin)
173 {
174  G4double tR = 0.895*CLHEP::fermi;
175  G4double pR = ParticleRadius(theParticle);
176 
177  G4double pZ = theParticle->GetPDGCharge()*fInvep;
178  G4double tZ = nucleon->GetPDGCharge()*fInvep;
179 
180  G4double pM = theParticle->GetPDGMass();
181  G4double tM = nucleon->GetPDGMass();
182 
183  G4double pElab = ekin + pM;
184  G4double totTcm = std::sqrt(pM*pM + tM*tM + 2.*pElab*tM) - pM -tM;
185 
186  G4double bC = fAlpha*pZ*tZ/(pR + tR);
187  return (totTcm > bC) ? 1. - bC/totTcm : 0.0;
188 }
189 
190 G4double G4NuclearRadii::CoulombFactor(
191  G4int Z, G4int A,
192  const G4ParticleDefinition* theParticle,
193  G4double ekin)
194 {
195  G4double tR = RadiusCB(Z, A);
196  G4double pR = ParticleRadius(theParticle);
197 
198  G4double pZ = theParticle->GetPDGCharge()*fInvep;
199 
200  G4double pM = theParticle->GetPDGMass();
201  G4double tM = G4NucleiProperties::GetNuclearMass(A, Z);
202 
203  G4double pElab = ekin + pM;
204  G4double totTcm = std::sqrt(pM*pM + tM*tM + 2.*pElab*tM) - pM -tM;
205 
206  G4double bC = fAlpha*pZ*Z/(pR + tR);
207  return (totTcm > bC) ? 1. - bC/totTcm : 0.0;
208 }
209 
210 const G4double G4NuclearRadii::r0[] = {
211  1.2,
212  1.3, 1.3, 1.3, 1.3,1.17,1.54,1.65,1.71, 1.7,1.75, // 1-10
213  1.7,1.57,1.53, 1.4, 1.3,1.30,1.44, 1.4, 1.4, 1.4, //11-20
214  1.4, 1.4,1.46, 1.4, 1.4,1.46,1.55, 1.5,1.38,1.48, //21-30
215  1.4, 1.4, 1.4,1.46, 1.4, 1.4, 1.4, 1.4, 1.4,1.45, //31-40
216  1.4, 1.4, 1.4, 1.4, 1.4, 1.4,1.45,1.48, 1.4,1.52, //41-50
217 1.46, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.5, //51-60
218  1.4, 1.4, 1.4, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 1.4, //61-70
219  1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3,1.33,1.43, //71-80
220  1.3,1.32,1.34, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, //81-90
221  1.3, 1.3};