ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GEMProbabilityVI.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4GEMProbabilityVI.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 // GEM de-excitation model
27 // by V. Ivanchenko (July 2019)
28 //
29 #include "G4GEMProbabilityVI.hh"
30 #include "G4NuclearLevelData.hh"
31 #include "G4LevelManager.hh"
32 #include "G4PairingCorrection.hh"
33 #include "G4NucleiProperties.hh"
34 #include "G4RandomDirection.hh"
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "Randomize.hh"
38 #include "G4Pow.hh"
39 #include "G4Exp.hh"
40 
41 // 10-Points Gauss-Legendre abcisas and weights
42 /*
43 const G4double G4GEMChannelVI::ws[] = {
44  0.0666713443086881,
45  0.149451349150581,
46  0.219086362515982,
47  0.269266719309996,
48  0.295524224714753,
49  0.295524224714753,
50  0.269266719309996,
51  0.219086362515982,
52  0.149451349150581,
53  0.0666713443086881
54  };
55 const G4double G4GEMChannelVI::xs[] = {
56  -0.973906528517172,
57  -0.865063366688985,
58  -0.679409568299024,
59  -0.433395394129247,
60  -0.148874338981631,
61  0.148874338981631,
62  0.433395394129247,
63  0.679409568299024,
64  0.865063366688985,
65  0.973906528517172
66 };
67 */
68 
69 G4GEMProbabilityVI::G4GEMProbabilityVI(G4int anA, G4int aZ, const G4LevelManager* p)
70  : G4VEmissionProbability(aZ, anA), lManager(p)
71 {
72  fragA = fragZ = 0;
73  resA13 = U = delta0 = delta1 = a0 = a1 = probmax = alphaP = betaP = 0.0;
74  Umax = bCoulomb = 0.0;
75  Gamma = 1.0;
76  pcoeff = Gamma*pEvapMass*CLHEP::millibarn
77  /((CLHEP::pi*CLHEP::hbarc)*(CLHEP::pi*CLHEP::hbarc));
78  coeff = CLHEP::fermi*CLHEP::fermi/(CLHEP::pi*CLHEP::hbarc*CLHEP::hbarc);
79 
80  isExcited = (!lManager || 0.0 == lManager->MaxLevelEnergy()) ? false : true;
81  A13 = pG4pow->Z13(theA);
82 
83  if(0 == aZ) {
84  ResetIntegrator(30, 0.25*CLHEP::MeV, 0.02);
85  } else {
86  ResetIntegrator(30, 0.5*CLHEP::MeV, 0.03);
87  }
88 }
89 
90 G4GEMProbabilityVI::~G4GEMProbabilityVI()
91 {}
92 
93 G4double G4GEMProbabilityVI::ComputeTotalProbability(
94  const G4Fragment& fragment, G4double CB)
95 {
96  fragA = fragment.GetA_asInt();
97  fragZ = fragment.GetZ_asInt();
98 
99  bCoulomb = CB;
100  U = fragment.GetExcitationEnergy();
101  delta0 = pNuclearLevelData->GetPairingCorrection(fragZ,fragA);
102  delta1 = pNuclearLevelData->GetPairingCorrection(resZ,resA);
103  Umax = pMass - pEvapMass - pResMass - CB;
104  if(0.0 >= Umax) { return 0.0; }
105 
106  resA13 = pG4pow->Z13(resA);
107  a0 = pNuclearLevelData->GetLevelDensity(fragZ,fragA,U);
108 
109  G4double C = 0.0;
110  G4int Z2 = theZ*theZ;
111  G4int Z3 = Z2*theZ;
112  G4int Z4 = Z2*Z2;
113 
114  if(resA >= 50) {
115  C = -0.10/(G4double)theA;
116  } else if(resZ > 20) {
117  C = (0.123482-0.00534691*theZ-0.0000610624*Z2+5.93719*1e-7*Z3+
118  1.95687*1e-8*Z4)/(G4double)theA;
119  }
120  if(0 == theZ) {
121  alphaP = 0.76+1.93/resA13;
122  betaP = (1.66/(resA13*resA13)-0.05)*CLHEP::MeV/alphaP;
123  } else {
124  alphaP = 1.0 + C;
125  betaP = - bCoulomb;
126  }
127  if(isExcited) {
128  pProbability = Integrated2DProbability();
129 
130  } else {
131  const G4double twoMass = pMass + pMass;
132  const G4double evapMass2 = pEvapMass*pEvapMass;
133  G4double ekinmax =
134  ((pMass-pResMass)*(pMass+pResMass) + evapMass2)/twoMass - pEvapMass;
135  G4double ekinmin =
136  std::max((CB*(twoMass - CB) + evapMass2)/twoMass - pEvapMass,0.0);
137  if(ekinmax <= ekinmin) { return 0.0; }
138  pProbability = IntegrateProbability(ekinmin, ekinmax, CB);
139  }
140  /*
141  G4cout << "G4GEMProbabilityVI: Z= " << theZ << " A= " << theA
142  << " resZ= " << resZ << " resA= " << resA
143  << " fragZ= " << fragZ << " fragA= " << fragA
144  << " prob= " << pProbability
145  << "\n U= " << U << " Umax= " << Umax << " d0= " << delta0
146  << " a0= " << a0 << G4endl;
147  */
148  return pProbability;
149 }
150 
151 G4double G4GEMProbabilityVI::ComputeProbability(G4double ekin, G4double)
152 {
153  // abnormal case - should never happens
154  if(pMass < pEvapMass + pResMass) { return 0.0; }
155 
156  const G4double m02 = pMass*pMass;
157  const G4double m12 = pEvapMass*pEvapMass;
158  const G4double mres = std::sqrt(m02 + m12 - 2.*pMass*(pEvapMass + ekin));
159 
160  G4double excRes = std::max(mres - pResMass, 0.0);
161  a1 = pNuclearLevelData->GetLevelDensity(resZ,resA,excRes);
162  G4double prob = ProbabilityDistributionFunction(0.0, excRes);
163 
164  //G4cout<<"### G4GEMProbabilityVI::ComputeProbability: Ekin(MeV)= "<<ekin
165  //<< " excRes(MeV)= " << excRes << " prob= " << prob << << G4endl;
166  return prob;
167 }
168 
169 G4Fragment* G4GEMProbabilityVI::SampleEvaporationFragment()
170 {
171  if(isExcited) { return Sample2DDistribution(); }
172  G4double ekin = SampleEnergy();
173  G4LorentzVector lv(std::sqrt(ekin*(ekin + 2.0*pEvapMass))
174  *G4RandomDirection(), ekin + pEvapMass);
175  G4Fragment* evFragment = new G4Fragment(theA, theZ, lv);
176  return evFragment;
177 }
178 
179 G4double G4GEMProbabilityVI::Integrated2DProbability()
180 {
181  return 0.0;
182 }
183 
184 G4double G4GEMProbabilityVI::ProbabilityDistributionFunction(
185  G4double exc, G4double resExc)
186 {
187  G4double Ux = (2.5 + 150.0/G4double(resA))*CLHEP::MeV;
188  G4double Ex = Ux + delta1;
189  G4double T = 1.0/(std::sqrt(a0/Ux) - 1.5/Ux);
190  G4double E0 = Ex - T*(G4Log(T) - G4Log(a0)*0.25
191  - 1.25*G4Log(Ux) + 2.0*std::sqrt(a0*Ux));
192 
193  G4double UxCN = (2.5 + 150.0/(G4double)theA)*CLHEP::MeV;
194  G4double ExCN = UxCN + delta0;
195  G4double TCN = 1.0/(std::sqrt(a0/UxCN) - 1.5/UxCN);
196 
197  G4double mass1 = pEvapMass + exc;
198  G4double mass2 = pResMass + resExc;
199 
200  G4double maxKinEnergy = std::max(0.5*((pMass - mass2)*(pMass + mass2)
201  + mass1*mass1)/pMass - mass1, 0.0);
202 
203  G4double Width = 0.0;
204  G4double t = maxKinEnergy/T;
205  if ( maxKinEnergy < Ex ) {
206  Width = (I1(t,t)*T + (betaP+bCoulomb)*I0(t))/G4Exp(E0/T);
207 
208  } else {
209 
210  G4double tx = Ex/T;
211  G4double s0 = 2.0*std::sqrt(a0*(maxKinEnergy-delta0));
212  G4double sx = 2.0*std::sqrt(a0*(Ex-delta0));
213 
214  // VI: protection against FPE exception
215  s0 = std::min(s0, 350.);
216 
217  G4double expE0T = G4Exp(E0/T);
218  G4double exps0 = G4Exp(s0);
219  const G4double sqrt2 = std::sqrt(2.0);
220 
221  Width = I1(t,tx)*T/expE0T + I3(s0,sx)*exps0/(sqrt2*a0);
222 
223  if (0 == theZ) {
224  Width += (betaP+bCoulomb)*(I0(tx)/expE0T + 2.0*sqrt2*I2(s0,sx)*exps0);
225  }
226  }
227  Width *= alphaP*pMass;
228 
229  //JMQ 190709 fix on Rb and geometrical cross sections according to
230  // Furihata's paper (JAERI-Data/Code 2001-105, p6)
231  G4double Rb = 0.0;
232  if (theA > 4) {
233  Rb = 1.12*(resA13 + A13) - 0.86*((resA13 + A13)/(resA13*A13))+2.85;
234  } else if (theA > 1) {
235  Rb=1.5*(resA13 + A13);
236  } else {
237  Rb = 1.5*resA13;
238  }
239 
240  G4double ild;
241  if (exc < ExCN ) {
242  G4double E0CN = ExCN - TCN*(G4Log(TCN) - 0.25*G4Log(a0)
243  - 1.25*G4Log(UxCN)
244  + 2.0*std::sqrt(a0*UxCN));
245  ild = G4Exp((exc-E0CN)/TCN)/TCN;
246  } else {
247  G4double x = exc - delta0;
248  G4double x1 = std::sqrt(a0*x);
249  ild = G4Exp(2*x1)/(x*std::sqrt(x1));
250  }
251 
252  Width *= (Rb*Rb/ild);
253  return Width;
254 }
255 
256 G4Fragment* G4GEMProbabilityVI::Sample2DDistribution()
257 {
258  G4Fragment* aFragment = nullptr;
259  return aFragment;
260 }
261 
262 G4double G4GEMProbabilityVI::I0(G4double t)
263 {
264  return G4Exp(t) - 1.0;
265 }
266 
267 G4double G4GEMProbabilityVI::I1(G4double t, G4double tx)
268 {
269  return (t - tx + 1.0)*G4Exp(tx) - t - 1.0;
270 }
271 
272 G4double G4GEMProbabilityVI::I2(G4double s0, G4double sx)
273 {
274  G4double S = 1.0/std::sqrt(s0);
275  G4double Sx = 1.0/std::sqrt(sx);
276 
277  G4double p1 = S*S*S*( 1.0 + S*S*( 1.5 + 3.75*S*S) );
278  G4double p2 = Sx*Sx*Sx*( 1.0 + Sx*Sx*( 1.5 + 3.75*Sx*Sx) )*G4Exp(sx-s0);
279 
280  return p1-p2;
281 }
282 
283 G4double G4GEMProbabilityVI::I3(G4double s0, G4double sx)
284 {
285  G4double s2 = s0*s0;
286  G4double sx2 = sx*sx;
287  G4double S = 1.0/std::sqrt(s0);
288  G4double S2 = S*S;
289  G4double Sx = 1.0/std::sqrt(sx);
290  G4double Sx2 = Sx*Sx;
291 
292  G4double p1 = S *(2.0 + S2 *( 4.0 + S2 *( 13.5 + S2 *( 60.0 + S2 * 325.125 ))));
293  G4double p2 = Sx*Sx2 *((s2-sx2) + Sx2 *((1.5*s2+0.5*sx2)
294  + Sx2 *((3.75*s2+0.25*sx2) + Sx2 *((12.875*s2+0.625*sx2)
295  + Sx2 *((59.0625*s2+0.9375*sx2) + Sx2 *(324.8*s2+3.28*sx2))))));
296  p2 *= G4Exp(sx-s0);
297  return p1-p2;
298 }
299