Geant4  10.02.p01
G4EvaporationProbability.cc
Go to the documentation of this file.
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 //J.M. Quesada (August2008). Based on:
27 //
28 // Hadronic Process: Nuclear De-excitations
29 // by V. Lara (Oct 1998)
30 //
31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse
32 // cross section option
33 // JMQ (06 September 2008) Also external choices have been added for
34 // superimposed Coulomb barrier (if useSICB is set true, by default is false)
35 //
36 // JMQ (14 february 2009) bug fixed in emission width: hbarc instead of
37 // hbar_Planck in the denominator
38 //
40 #include "G4VCoulombBarrier.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4PairingCorrection.hh"
44 #include "G4NucleiProperties.hh"
45 #include "G4KalbachCrossSection.hh"
47 #include "Randomize.hh"
48 #include "G4Exp.hh"
49 
50 using namespace std;
51 
52 static const G4double explim = 350.;
53 static const G4double invmev = 1./CLHEP::MeV;
54 static const G4double ssqr3 = 1.5*std::sqrt(3.0);
55 
57  G4double aGamma,
59  : theA(anA),
60  theZ(aZ),
61  Gamma(aGamma)
62 {
63  resZ = resA = fragA = fragZ = nbins = 0;
64  resA13 = muu = resMass = fragMass = U = delta0 = delta1 = a0 = 0.0;
66  index = 0;
67  if(1 == theZ) { index = theA; }
68  else { index = theA + 1; }
69  for(G4int i=0; i<11; ++i) { probability[i] = 0.0; }
70 }
71 
73 {}
74 
75 // obsolete method
77  const G4Fragment&, G4double)
78 {
79  return 0.0;
80 }
81 
83  const G4Fragment & fragment, G4double minEnergy, G4double maxEnergy)
84 {
85  if (maxEnergy <= minEnergy) { return 0.0; }
86 
87  fragA = fragment.GetA_asInt();
88  fragZ = fragment.GetZ_asInt();
89  resA = fragA - theA;
90  resZ = fragZ - theZ;
91 
92  U = fragment.GetExcitationEnergy();
94  if(U < delta0) { return 0.0; }
95 
96  fragMass = fragment.GetGroundStateMass();
99  resA13 = fG4pow->Z13(resA);
100 
101  G4double Width = 0.0;
102  if (OPTxs==0) {
103 
104  G4double SystemEntropy = 2.0*std::sqrt(
106 
107  static const G4double RN2 =
108  2.25*fermi*fermi/(twopi* hbar_Planck*hbar_Planck);
109 
110  G4double Alpha = CalcAlphaParam(fragment);
111  G4double Beta = CalcBetaParam(fragment);
112 
114  G4double GlobalFactor = Gamma*Alpha*partMass*RN2*resA13*resA13/(a0*a0);
115 
116  G4double maxea = maxEnergy*a0;
117  G4double Term1 = Beta*a0 - 1.5 + maxea;
118  G4double Term2 = (2.0*Beta*a0-3.0)*std::sqrt(maxea) + 2*maxea;
119 
120  G4double ExpTerm1 = 0.0;
121  if (SystemEntropy <= explim) { ExpTerm1 = G4Exp(-SystemEntropy); }
122 
123  G4double ExpTerm2 = 2.*std::sqrt(maxea) - SystemEntropy;
124  ExpTerm2 = std::min(ExpTerm2, explim);
125  ExpTerm2 = G4Exp(ExpTerm2);
126 
127  Width = GlobalFactor*(Term1*ExpTerm1 + Term2*ExpTerm2);
128 
129  } else {
130 
132  // compute power once
133  if(OPTxs <= 2) {
135  } else {
137  }
138  // if Coulomb barrier cutoff is superimposed for all cross sections
139  // then the limit is the Coulomb Barrier
140  Width = IntegrateEmissionProbability(minEnergy, maxEnergy);
141  }
142  return Width;
143 }
144 
145 G4double
147 {
148  G4double sum = 0.0;
149  G4double del = up - low;
150  nbins = std::min(10, G4int(del*invmev));
151  nbins = std::max(nbins, 2);
152  del /= G4double(nbins);
153  G4double T = low - del*0.5;
154  for(G4int i=1; i<=nbins; ++i) {
155  T += del;
157  probability[i] = sum;
158  }
159  sum *= del;
160  return sum;
161 }
162 
164 {
165  //G4cout << "### G4EvaporationProbability::ProbabilityDistributionFunction"
166  // << G4endl;
167 
168  G4double E0 = U - delta0;
169  G4double E1 = U + fragMass - partMass - resMass - delta1 - K;
170  /*
171  G4cout << "PDF: FragZ= " << fragment.GetZ_asInt() << " FragA= "
172  << fragment.GetA_asInt()
173  << " Z= " << theZ << " A= " << theA
174  << " K= " << K << " E0= " << E0 << " E1= " << E1 << G4endl;
175  */
176  if(E1 < 0.0) { return 0.0; }
177 
179 
180  //JMQ 14/02/09 BUG fixed: hbarc should be in the denominator instead
181  // of hbar_Planck; without 1/hbar_Panck remains as a width
182 
183  static const G4double pcoeff = millibarn/((pi*hbarc)*(pi*hbarc));
184 
185  // Fixed numerical problem
186  G4double Prob = pcoeff*Gamma*partMass
187  *G4Exp(2*(std::sqrt(a1*E1) - std::sqrt(a0*E0)))*K*CrossSection(K);
188 
189  return Prob;
190 }
191 
192 G4double
194 {
195  G4double res;
196  if(OPTxs <= 2) {
198  index, theZ,
199  resZ, resA);
200  } else {
202  index, theZ, theA,
203  resZ, resA);
204  }
205  return res;
206 }
207 
208 G4double
210  G4double maxKinEnergy)
211 {
212  if(maxKinEnergy <= minKinEnergy) { return 0.0; }
213 
214  G4double T = 0.0;
215  if (OPTxs==0) {
216  // JMQ:
217  // It uses Dostrovsky's approximation for the inverse reaction cross
218  // in the probability for fragment emission
219  // MaximalKineticEnergy energy in the original version (V.Lara) was
220  // calculated at the Coulomb barrier.
221 
222  G4double Rb = 4.0*a0*maxKinEnergy;
223  G4double RbSqrt = std::sqrt(Rb);
224  G4double PEX1 = 0.0;
225  if (RbSqrt < 160.0) PEX1 = G4Exp(-RbSqrt);
226  G4double Rk = 0.0;
227  G4double FRk = 0.0;
228  do {
229  G4double RandNumber = G4UniformRand();
230  Rk = 1.0 + (1./RbSqrt)*G4Log(RandNumber + (1.0-RandNumber)*PEX1);
231  G4double Q1 = 1.0;
232  G4double Q2 = 1.0;
233  if (theZ == 0) { // for emitted neutron
234  G4double Beta = (2.12/(resA*resA) - 0.05)*MeV/(0.76 + 2.2/resA13);
235  Q1 = 1.0 + Beta/maxKinEnergy;
236  Q2 = Q1*std::sqrt(Q1);
237  }
238 
239  FRk = ssqr3 * Rk * (Q1 - Rk*Rk)/Q2;
240 
241  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
242  } while (FRk < G4UniformRand());
243 
244  T = maxKinEnergy * (1.0-Rk*Rk) + minKinEnergy;
245 
246  } else {
247 
249  G4int i = 0;
250  for(; i<nbins; ++i) { if(p <= probability[i+1]) { break; } }
251  G4double delta = (maxKinEnergy - minKinEnergy)/G4double(nbins);
252  T = minKinEnergy + delta*i
253  + delta*(p - probability[i])/(probability[i+1] - probability[i]);
254 
255  }
256  //G4cout << " T= " << T << G4endl;
257  return T;
258 }
static const double MeV
Definition: G4SIunits.hh:211
static G4double GetNuclearMass(const G4double A, const G4double Z)
static const G4double ssqr3
static const G4double invmev
static const G4double a1
static G4double ComputePowerParameter(G4int resA, G4int idx)
G4double SampleKineticEnergy(G4double minKineticEnergy, G4double maxKineticEnergy)
static G4double ComputePowerParameter(G4int resA, G4int idx)
static G4double ComputeCrossSection(G4double K, G4double resA13, G4double amu1, G4int idx, G4int Z, G4int resZ, G4int resA)
int G4int
Definition: G4Types.hh:78
G4EvaporationLevelDensityParameter * theEvapLDPptr
#define G4UniformRand()
Definition: Randomize.hh:97
G4double Z13(G4int Z) const
Definition: G4Pow.hh:127
G4int GetA_asInt() const
Definition: G4Fragment.hh:251
static const G4double explim
static const double twopi
Definition: G4SIunits.hh:75
virtual G4double CalcBetaParam(const G4Fragment &fragment)=0
G4double GetPairingCorrection(G4int A, G4int Z) const
G4double EmissionProbability(const G4Fragment &fragment, G4double maxKineticEnergy)
G4double LevelDensityParameter(G4int A, G4int Z, G4double U) const
G4EvaporationProbability(G4int anA, G4int aZ, G4double aGamma, G4VCoulombBarrier *)
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:273
virtual G4double CalcAlphaParam(const G4Fragment &fragment)=0
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static const double pi
Definition: G4SIunits.hh:74
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static const double millibarn
Definition: G4SIunits.hh:105
G4int GetZ_asInt() const
Definition: G4Fragment.hh:256
G4PairingCorrection * fPairCorr
double G4double
Definition: G4Types.hh:76
G4double IntegrateEmissionProbability(G4double low, G4double up)
static G4double ComputeCrossSection(G4double K, G4double resA13, G4double amu1, G4int idx, G4int Z, G4int A, G4int resZ, G4int resA)
G4double ProbabilityDistributionFunction(G4double K)
G4double TotalProbability(const G4Fragment &fragment, G4double minKineticEnergy, G4double maxKineticEnergy)
static const double fermi
Definition: G4SIunits.hh:102
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:268