Geant4  10.01.p02
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 hbar_Planck in the denominator
37 //
38 #include <iostream>
39 
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4PairingCorrection.hh"
44 #include "G4ParticleTable.hh"
45 #include "G4IonTable.hh"
46 
47 using namespace std;
48 
50  G4double aGamma,
51  G4VCoulombBarrier * aCoulombBarrier)
52  : theA(anA),
53  theZ(aZ),
54  Gamma(aGamma),
55  theCoulombBarrierptr(aCoulombBarrier)
56 {}
57 
59  : theA(0),
60  theZ(0),
61  Gamma(0.0),
62  theCoulombBarrierptr(0)
63 {}
64 
66 {}
67 
68 G4double
70  G4double anEnergy)
71 {
72  G4double probability = 0.0;
73 
74  if (anEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) {
75  probability = CalculateProbability(fragment, anEnergy);
76 
77  }
78  return probability;
79 }
80 
82 
83 // Computes the integrated probability of evaporation channel
84 G4double
86  G4double MaximalKineticEnergy)
87 {
88  G4int ResidualA = fragment.GetA_asInt() - theA;
89  G4int ResidualZ = fragment.GetZ_asInt() - theZ;
90  G4double U = fragment.GetExcitationEnergy();
91 
92  if (OPTxs==0) {
93 
94  G4double NuclearMass = fragment.ComputeGroundStateMass(theZ,theA);
95 
97  fragment.GetZ_asInt());
98 
99  G4double SystemEntropy = 2.0*std::sqrt(
100  theEvapLDPptr->LevelDensityParameter(fragment.GetA_asInt(),fragment.GetZ_asInt(),U)*
101  (U-delta0));
102 
103  static const G4double RN = 1.5*fermi;
104 
105  G4double Alpha = CalcAlphaParam(fragment);
106  G4double Beta = CalcBetaParam(fragment);
107 
108  G4double Rmax = MaximalKineticEnergy;
109  G4double a = theEvapLDPptr->LevelDensityParameter(ResidualA,ResidualZ,Rmax);
110  G4double GlobalFactor = Gamma * Alpha/(a*a) *
111  (NuclearMass*RN*RN*fG4pow->Z23(ResidualA))/
112  (twopi* hbar_Planck*hbar_Planck);
113  G4double Term1 = (2.0*Beta*a-3.0)/2.0 + Rmax*a;
114  G4double Term2 = (2.0*Beta*a-3.0)*std::sqrt(Rmax*a) + 2.0*a*Rmax;
115 
116  G4double ExpTerm1 = 0.0;
117  if (SystemEntropy <= 600.0) { ExpTerm1 = std::exp(-SystemEntropy); }
118 
119  G4double ExpTerm2 = 2.*std::sqrt(a*Rmax) - SystemEntropy;
120  if (ExpTerm2 > 700.0) { ExpTerm2 = 700.0; }
121  ExpTerm2 = std::exp(ExpTerm2);
122 
123  G4double Width = GlobalFactor*(Term1*ExpTerm1 + Term2*ExpTerm2);
124 
125  return Width;
126 
127  } else if (OPTxs==1 || OPTxs==2 ||OPTxs==3 || OPTxs==4) {
128 
129  G4double EvaporatedMass = fragment.ComputeGroundStateMass(theZ,theA);
130  G4double ResidulalMass = fragment.ComputeGroundStateMass(ResidualZ,ResidualA);
131  G4double limit = std::max(0.0,fragment.GetGroundStateMass()-EvaporatedMass-ResidulalMass);
132  if (useSICB) {
133  limit = std::max(limit,theCoulombBarrierptr->GetCoulombBarrier(ResidualA,ResidualZ,U));
134  }
135 
136  if (MaximalKineticEnergy <= limit) { return 0.0; }
137 
138  // if Coulomb barrier cutoff is superimposed for all cross sections
139  // then the limit is the Coulomb Barrier
140  G4double LowerLimit= limit;
141 
142  //MaximalKineticEnergy: asimptotic value (already accounted for in G4EvaporationChannel)
143 
144  G4double UpperLimit = MaximalKineticEnergy;
145 
146  G4double Width = IntegrateEmissionProbability(fragment,LowerLimit,UpperLimit);
147 
148  return Width;
149  } else {
150  std::ostringstream errOs;
151  errOs << "Bad option for cross sections at evaporation" <<G4endl;
152  throw G4HadronicException(__FILE__, __LINE__, errOs.str());
153  }
154 
155 }
156 
158 
161  const G4double & Low, const G4double & Up )
162 {
163  static const G4int N = 10;
164  // 10-Points Gauss-Legendre abcisas and weights
165  static const G4double w[N] = {
166  0.0666713443086881,
167  0.149451349150581,
168  0.219086362515982,
169  0.269266719309996,
170  0.295524224714753,
171  0.295524224714753,
172  0.269266719309996,
173  0.219086362515982,
174  0.149451349150581,
175  0.0666713443086881
176  };
177  static const G4double x[N] = {
178  -0.973906528517172,
179  -0.865063366688985,
180  -0.679409568299024,
181  -0.433395394129247,
182  -0.148874338981631,
183  0.148874338981631,
184  0.433395394129247,
185  0.679409568299024,
186  0.865063366688985,
187  0.973906528517172
188  };
189 
190  G4double Total = 0.0;
191 
192 
193  for (G4int i = 0; i < N; i++)
194  {
195 
196  G4double KineticE = ((Up-Low)*x[i]+(Up+Low))/2.0;
197 
198  Total += w[i]*ProbabilityDistributionFunction(fragment, KineticE);
199 
200  }
201  Total *= (Up-Low)/2.0;
202  return Total;
203 }
204 
205 
207 //New method (OPT=1,2,3,4)
208 
209 G4double
211  G4double K)
212 {
213  G4int ResidualA = fragment.GetA_asInt() - theA;
214  G4int ResidualZ = fragment.GetZ_asInt() - theZ;
215  G4double U = fragment.GetExcitationEnergy();
216  //G4cout << "### G4EvaporationProbability::ProbabilityDistributionFunction" << G4endl;
217  //G4cout << "FragZ= " << fragment.GetZ_asInt() << " FragA= " << fragment.GetA_asInt()
218  // << " Z= " << theZ << " A= " << theA << G4endl;
219  //G4cout << "PC " << fPairCorr << " DP " << theEvapLDPptr << G4endl;
220 
221  // if(K <= theCoulombBarrierptr->GetCoulombBarrier(ResidualA,ResidualZ,U)) return 0.0;
222 
223  G4double delta1 = fPairCorr->GetPairingCorrection(ResidualA,ResidualZ);
224 
225  G4double delta0 = fPairCorr->GetPairingCorrection(fragment.GetA_asInt(),
226  fragment.GetZ_asInt());
227 
228 
229  G4double ParticleMass = fragment.ComputeGroundStateMass(theZ,theA);
230  G4double ResidualMass = fragment.ComputeGroundStateMass(ResidualZ,ResidualA);
231 
232  G4double theSeparationEnergy = ParticleMass + ResidualMass
233  - fragment.GetGroundStateMass();
234 
236  fragment.GetZ_asInt(),
237  U - delta0);
238 
239  G4double a1 = theEvapLDPptr->LevelDensityParameter(ResidualA, ResidualZ,
240  U - theSeparationEnergy - delta1);
241 
242 
243  G4double E0 = U - delta0;
244 
245  G4double E1 = U - theSeparationEnergy - delta1 - K;
246 
247  if (E1<0.) { return 0.; }
248 
249  //JMQ 14/02/09 BUG fixed: hbarc should be in the denominator instead of hbar_Planck
250  //Without 1/hbar_Panck remains as a width
251 
252  //G4double Prob=Gamma*ParticleMass/((pi*hbarc)*(pi*hbarc)*std::exp(2*std::sqrt(a0*E0)))
253  // *K*CrossSection(fragment,K)*std::exp(2*std::sqrt(a1*E1))*millibarn;
254 
255  static const G4double pcoeff = millibarn/((pi*hbarc)*(pi*hbarc));
256 
257  // Fixed numerical problem
258  G4double Prob = pcoeff*Gamma*ParticleMass*std::exp(2*(std::sqrt(a1*E1) - std::sqrt(a0*E0)))
259  *K*CrossSection(fragment,K);
260 
261  return Prob;
262 }
263 
264 
const G4double a0
G4double EmissionProbability(const G4Fragment &fragment, G4double anEnergy)
static const G4double a1
const G4double pi
G4double a
Definition: TRTMaterials.hh:39
G4double IntegrateEmissionProbability(const G4Fragment &aFragment, const G4double &Low, const G4double &Up)
int G4int
Definition: G4Types.hh:78
G4double ProbabilityDistributionFunction(const G4Fragment &aFragment, G4double K)
G4EvaporationLevelDensityParameter * theEvapLDPptr
G4int GetA_asInt() const
Definition: G4Fragment.hh:243
virtual G4double CalcBetaParam(const G4Fragment &fragment)=0
G4double GetPairingCorrection(G4int A, G4int Z) const
G4double LevelDensityParameter(G4int A, G4int Z, G4double U) const
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:265
virtual G4double CalcAlphaParam(const G4Fragment &fragment)=0
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static const double millibarn
Definition: G4SIunits.hh:96
G4int GetZ_asInt() const
Definition: G4Fragment.hh:248
G4double Z23(G4int Z) const
Definition: G4Pow.hh:151
virtual G4double GetCoulombBarrier(G4int ARes, G4int ZRes, G4double U) const =0
#define G4endl
Definition: G4ios.hh:61
G4PairingCorrection * fPairCorr
virtual G4double CrossSection(const G4Fragment &fragment, G4double K)=0
double G4double
Definition: G4Types.hh:76
G4VCoulombBarrier * theCoulombBarrierptr
G4double ComputeGroundStateMass(G4int Z, G4int A) const
Definition: G4Fragment.hh:298
static const double fermi
Definition: G4SIunits.hh:93
G4double CalculateProbability(const G4Fragment &fragment, G4double MaximalKineticEnergy)
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:260