| Geant4 Cross Reference |
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 G4GEMProbability
30 //
31 //
32 // Hadronic Process: Nuclear De-excitations
33 // by V. Lara (Sept 2001)
34 //
35 //
36 // Hadronic Process: Nuclear De-excitations
37 // by V. Lara (Sept 2001)
38 //
39 // J. M. Quesada : several fixes in total GEM width
40 // J. M. Quesada 14/07/2009 bug fixed in total emission width (hbarc)
41 // J. M. Quesada (September 2009) several fixes:
42 // -level density parameter of residual calculated at its right excitation energy.
43 // -InitialLeveldensity calculated according to the right conditions of the
44 // initial excited nucleus.
45 // J. M. Quesada 19/04/2010 fix in emission probability calculation.
46 // V.Ivanchenko 20/04/2010 added usage of G4Pow and use more safe computation
47 // V.Ivanchenko 18/05/2010 trying to speedup the most slow method
48 // by usage of G4Pow, integer Z and A; moved constructor,
49 // destructor and virtual functions to source
50
51 #include "G4GEMProbability.hh"
52 #include "G4PairingCorrection.hh"
53 #include "G4NucleiProperties.hh"
54 #include "G4PhysicalConstants.hh"
55 #include "G4SystemOfUnits.hh"
56 #include "G4Log.hh"
57
58 G4GEMProbability:: G4GEMProbability(G4int anA, G4int aZ, G4double aSpin)
59 : G4VEmissionProbability(aZ, anA), Spin(aSpin),
60 theCoulombBarrierPtr(nullptr)
61 {
62 theEvapLDPptr = new G4EvaporationLevelDensityParameter;
63 fG4pow = G4Pow::GetInstance();
64 fPlanck= CLHEP::hbar_Planck*fG4pow->logZ(2);
65 fNucData = G4NuclearLevelData::GetInstance();
66 }
67
68 G4GEMProbability::~G4GEMProbability()
69 {
70 delete theEvapLDPptr;
71 }
72
73 G4double G4GEMProbability::EmissionProbability(const G4Fragment & fragment,
74 G4double MaximalKineticEnergy)
75 {
76 G4double probability = 0.0;
77
78 if (MaximalKineticEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) {
79 G4double CoulombBarrier = GetCoulombBarrier(fragment);
80 G4double InitialLevelDensity = ComputeInitialLevelDensity(fragment);
81 G4double Ux, UxSqrt, UxLog;
82 PrecomputeResidualQuantities(fragment, Ux, UxSqrt, UxLog);
83
84 probability =
85 CalcProbability(fragment,MaximalKineticEnergy,CoulombBarrier,Spin,
86 InitialLevelDensity,Ux,UxSqrt,UxLog);
87
88 // Next there is a loop over excited states for this channel
89 // summing probabilities
90 std::size_t nn = ExcitEnergies.size();
91 if (0 < nn) {
92 for (std::size_t i = 0; i <nn; ++i) {
93 // substract excitation energies
94 G4double Tmax = MaximalKineticEnergy - ExcitEnergies[i];
95 if (Tmax > 0.0) {
96 G4double width = CalcProbability(fragment,Tmax,CoulombBarrier,ExcitSpins[i],
97 InitialLevelDensity,Ux,UxSqrt,UxLog);
98 //JMQ April 2010 added condition to prevent reported crash
99 // update probability
100 if (width > 0. && fPlanck < width*ExcitLifetimes[i]) {
101 probability += width;
102 }
103 }
104 }
105 }
106 }
107 // Normalization = probability;
108 return probability;
109 }
110
111 G4double G4GEMProbability::ComputeInitialLevelDensity(const G4Fragment & fragment) const
112 {
113 G4int A = fragment.GetA_asInt();
114 G4int Z = fragment.GetZ_asInt();
115 G4double U = fragment.GetExcitationEnergy();
116
117 // ***PARENT***
118 //JMQ (September 2009) the following quantities refer to the PARENT:
119
120 G4double deltaCN = fNucData->GetPairingCorrection(Z, A);
121 G4double aCN = theEvapLDPptr->LevelDensityParameter(A, Z, U-deltaCN);
122 G4double UxCN = (2.5 + 150.0/G4double(A))*MeV;
123 G4double ExCN = UxCN + deltaCN;
124 G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
125 // ***end PARENT***
126
127 //JMQ 160909 fix: initial level density must be calculated according to the
128 // conditions at the initial compound nucleus
129 // (it has been removed from previous "if" for the residual)
130
131 G4double InitialLevelDensity;
132 if ( U < ExCN )
133 {
134 //JMQ fixed bug in units
135 //VI moved the computation here
136 G4double E0CN = ExCN - TCN*(G4Log(TCN/MeV) - 0.25*G4Log(aCN*MeV)
137 - 1.25*G4Log(UxCN/MeV)
138 + 2.0*std::sqrt(aCN*UxCN));
139
140 InitialLevelDensity = (pi/12.0)*G4Exp((U-E0CN)/TCN)/TCN;
141 }
142 else
143 {
144 //VI speedup
145 G4double x = U-deltaCN;
146 G4double x1 = std::sqrt(aCN*x);
147
148 InitialLevelDensity = (pi/12.0)*G4Exp(2*x1)/(x*std::sqrt(x1));
149 }
150
151 return InitialLevelDensity;
152 }
153
154 void G4GEMProbability::PrecomputeResidualQuantities(const G4Fragment & fragment,
155 G4double &Ux,
156 G4double &UxSqrt,
157 G4double &UxLog) const
158 {
159 G4int A = fragment.GetA_asInt();
160 G4int ResidualA = A - theA;
161
162 Ux = (2.5 + 150.0/G4double(ResidualA))*MeV;
163 UxSqrt = std::sqrt(Ux);
164 UxLog = G4Log(Ux/MeV);
165 }
166
167 G4double G4GEMProbability::CalcProbability(const G4Fragment & fragment,
168 G4double MaximalKineticEnergy,
169 G4double V, G4double spin,
170 G4double InitialLevelDensity,
171 G4double Ux, G4double UxSqrt,
172 G4double UxLog) const
173
174 // Calculate integrated probability (width) for evaporation channel
175 {
176 G4int A = fragment.GetA_asInt();
177 G4int Z = fragment.GetZ_asInt();
178
179 G4int ResidualA = A - theA;
180 G4int ResidualZ = Z - theZ;
181
182 G4double NuclearMass = fragment.ComputeGroundStateMass(theZ, theA);
183
184 G4double Alpha = CalcAlphaParam(fragment);
185 G4double Beta = CalcBetaParam(fragment);
186
187 // ***RESIDUAL***
188 //JMQ (September 2009) the following quantities refer to the RESIDUAL:
189
190 G4double delta0 = fNucData->GetPairingCorrection(ResidualZ, ResidualA);
191
192 G4double a = theEvapLDPptr->
193 LevelDensityParameter(ResidualA,ResidualZ,MaximalKineticEnergy+V-delta0);
194 G4double aSqrt = std::sqrt(a);
195 G4double Ex = Ux + delta0;
196 G4double T = 1.0/(aSqrt/UxSqrt - 1.5/Ux);
197 //JMQ fixed bug in units
198 G4double E0 = Ex - T*(G4Log(T/MeV) - G4Log(a*MeV)/4.0
199 - 1.25*UxLog + 2.0*aSqrt*UxSqrt);
200 // ***end RESIDUAL ***
201
202 G4double Width;
203 G4double t = MaximalKineticEnergy/T;
204 if ( MaximalKineticEnergy < Ex ) {
205 //JMQ 190709 bug in I1 fixed (T was missing)
206 Width = (I1(t,t)*T + (Beta+V)*I0(t))/G4Exp(E0/T);
207 } else {
208
209 //VI minor speedup
210 G4double expE0T = G4Exp(E0/T);
211 static const G4double sqrt2 = std::sqrt(2.0);
212
213 G4double tx = Ex/T;
214 G4double s0 = 2.0*std::sqrt(a*(MaximalKineticEnergy-delta0));
215 G4double sx = 2.0*std::sqrt(a*(Ex-delta0));
216 // VI: protection against FPE exception
217 if(s0 > 350.) { s0 = 350.; }
218 Width = I1(t,tx)*T/expE0T + I3(s0,sx)*G4Exp(s0)/(sqrt2*a);
219
220 // VI this cannot happen!
221 // For charged particles (Beta+V) = 0 beacuse Beta = -V
222 //if (theZ == 0) {
223 // Width += (Beta+V)*(I0(tx)/expE0T + 2.0*sqrt2*I2(s0,sx)*G4Exp(s0));
224 //}
225 }
226
227 //JMQ 14/07/2009 BIG BUG : NuclearMass is in MeV => hbarc instead of
228 // hbar_planck must be used
229 // G4double g = (2.0*spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
230 G4double gg = (2.0*spin+1.0)*NuclearMass/(pi2* hbarc*hbarc);
231
232 //JMQ 190709 fix on Rb and geometrical cross sections according to
233 // Furihata's paper (JAERI-Data/Code 2001-105, p6)
234 G4double Rb = 0.0;
235 if (theA > 4)
236 {
237 G4double Ad = fG4pow->Z13(ResidualA);
238 G4double Aj = fG4pow->Z13(theA);
239 Rb = 1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85;
240 Rb *= fermi;
241 }
242 else if (theA>1)
243 {
244 G4double Ad = fG4pow->Z13(ResidualA);
245 G4double Aj = fG4pow->Z13(theA);
246 Rb=1.5*(Aj+Ad)*fermi;
247 }
248 else
249 {
250 G4double Ad = fG4pow->Z13(ResidualA);
251 Rb = 1.5*Ad*fermi;
252 }
253 G4double GeometricalXS = pi*Rb*Rb;
254 //end of JMQ fix on Rb by 190709
255
256 //JMQ 190709 BUG : pi instead of sqrt(pi) must be here according
257 // to Furihata's report:
258 Width *= pi*gg*GeometricalXS*Alpha/(12.0*InitialLevelDensity);
259
260 return Width;
261 }
262
263 G4double G4GEMProbability::I3(G4double s0, G4double sx) const
264 {
265 G4double s2 = s0*s0;
266 G4double sx2 = sx*sx;
267 G4double S = 1.0/std::sqrt(s0);
268 G4double S2 = S*S;
269 G4double Sx = 1.0/std::sqrt(sx);
270 G4double Sx2 = Sx*Sx;
271
272 G4double p1 = S *(2.0 + S2 *( 4.0 + S2 *( 13.5 + S2 *( 60.0 + S2 * 325.125 ))));
273 G4double p2 = Sx*Sx2 *(
274 (s2-sx2) + Sx2 *(
275 (1.5*s2+0.5*sx2) + Sx2 *(
276 (3.75*s2+0.25*sx2) + Sx2 *(
277 (12.875*s2+0.625*sx2) + Sx2 *(
278 (59.0625*s2+0.9375*sx2) + Sx2 *(324.8*s2+3.28*sx2))))));
279
280 p2 *= G4Exp(sx-s0);
281
282 return p1-p2;
283 }
284
285 void G4GEMProbability::Dump() const
286 {
287 G4double mass = G4NucleiProperties::GetNuclearMass(theA, theZ);
288 G4double efermi = 0.0;
289 if(theA > 1) {
290 efermi = G4NucleiProperties::GetNuclearMass(theA-1, theZ)
291 + neutron_mass_c2 - mass;
292 }
293 std::size_t nlev = ExcitEnergies.size();
294 G4cout << "GEM: List of Excited States for Isotope Z= "
295 << theZ << " A= " << theA << " Nlevels= " << nlev
296 << " Efermi(MeV)= " << efermi
297 << G4endl;
298 for(std::size_t i=0; i< nlev; ++i) {
299 G4cout << "Z= " << theZ << " A= " << theA
300 << " Mass(GeV)= " << mass/GeV
301 << " Eexc(MeV)= " << ExcitEnergies[i]
302 << " Time(ns)= " << ExcitLifetimes[i]/ns
303 << G4endl;
304 }
305 G4cout << G4endl;
306 }
307