00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038 #include <iostream>
00039
00040 #include "G4EvaporationProbability.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4PairingCorrection.hh"
00044 #include "G4ParticleTable.hh"
00045 #include "G4IonTable.hh"
00046
00047 using namespace std;
00048
00049 G4EvaporationProbability::G4EvaporationProbability(G4int anA, G4int aZ,
00050 G4double aGamma,
00051 G4VCoulombBarrier * aCoulombBarrier)
00052 : theA(anA),
00053 theZ(aZ),
00054 Gamma(aGamma),
00055 theCoulombBarrierptr(aCoulombBarrier)
00056 {}
00057
00058 G4EvaporationProbability::G4EvaporationProbability()
00059 : theA(0),
00060 theZ(0),
00061 Gamma(0.0),
00062 theCoulombBarrierptr(0)
00063 {}
00064
00065 G4EvaporationProbability::~G4EvaporationProbability()
00066 {}
00067
00068 G4double
00069 G4EvaporationProbability::EmissionProbability(const G4Fragment & fragment,
00070 G4double anEnergy)
00071 {
00072 G4double probability = 0.0;
00073
00074 if (anEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) {
00075 probability = CalculateProbability(fragment, anEnergy);
00076
00077 }
00078 return probability;
00079 }
00080
00082
00083
00084 G4double
00085 G4EvaporationProbability::CalculateProbability(const G4Fragment & fragment,
00086 G4double MaximalKineticEnergy)
00087 {
00088 G4int ResidualA = fragment.GetA_asInt() - theA;
00089 G4int ResidualZ = fragment.GetZ_asInt() - theZ;
00090 G4double U = fragment.GetExcitationEnergy();
00091
00092 if (OPTxs==0) {
00093
00094 G4double NuclearMass = fragment.ComputeGroundStateMass(theZ,theA);
00095
00096 G4double delta0 = fPairCorr->GetPairingCorrection(fragment.GetA_asInt(),
00097 fragment.GetZ_asInt());
00098
00099 G4double SystemEntropy = 2.0*std::sqrt(
00100 theEvapLDPptr->LevelDensityParameter(fragment.GetA_asInt(),fragment.GetZ_asInt(),U)*
00101 (U-delta0));
00102
00103 const G4double RN = 1.5*fermi;
00104
00105 G4double Alpha = CalcAlphaParam(fragment);
00106 G4double Beta = CalcBetaParam(fragment);
00107
00108 G4double Rmax = MaximalKineticEnergy;
00109 G4double a = theEvapLDPptr->LevelDensityParameter(ResidualA,ResidualZ,Rmax);
00110 G4double GlobalFactor = Gamma * Alpha/(a*a) *
00111 (NuclearMass*RN*RN*fG4pow->Z23(ResidualA))/
00112 (twopi* hbar_Planck*hbar_Planck);
00113 G4double Term1 = (2.0*Beta*a-3.0)/2.0 + Rmax*a;
00114 G4double Term2 = (2.0*Beta*a-3.0)*std::sqrt(Rmax*a) + 2.0*a*Rmax;
00115
00116 G4double ExpTerm1 = 0.0;
00117 if (SystemEntropy <= 600.0) { ExpTerm1 = std::exp(-SystemEntropy); }
00118
00119 G4double ExpTerm2 = 2.*std::sqrt(a*Rmax) - SystemEntropy;
00120 if (ExpTerm2 > 700.0) { ExpTerm2 = 700.0; }
00121 ExpTerm2 = std::exp(ExpTerm2);
00122
00123 G4double Width = GlobalFactor*(Term1*ExpTerm1 + Term2*ExpTerm2);
00124
00125 return Width;
00126
00127 } else if (OPTxs==1 || OPTxs==2 ||OPTxs==3 || OPTxs==4) {
00128
00129 G4double EvaporatedMass = fragment.ComputeGroundStateMass(theZ,theA);
00130 G4double ResidulalMass = fragment.ComputeGroundStateMass(ResidualZ,ResidualA);
00131 G4double limit = std::max(0.0,fragment.GetGroundStateMass()-EvaporatedMass-ResidulalMass);
00132 if (useSICB) {
00133 limit = std::max(limit,theCoulombBarrierptr->GetCoulombBarrier(ResidualA,ResidualZ,U));
00134 }
00135
00136 if (MaximalKineticEnergy <= limit) { return 0.0; }
00137
00138
00139
00140 G4double LowerLimit= limit;
00141
00142
00143
00144 G4double UpperLimit = MaximalKineticEnergy;
00145
00146 G4double Width = IntegrateEmissionProbability(fragment,LowerLimit,UpperLimit);
00147
00148 return Width;
00149 } else {
00150 std::ostringstream errOs;
00151 errOs << "Bad option for cross sections at evaporation" <<G4endl;
00152 throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00153 }
00154
00155 }
00156
00158
00159 G4double G4EvaporationProbability::
00160 IntegrateEmissionProbability(const G4Fragment & fragment,
00161 const G4double & Low, const G4double & Up )
00162 {
00163 static const G4int N = 10;
00164
00165 static const G4double w[N] = {
00166 0.0666713443086881,
00167 0.149451349150581,
00168 0.219086362515982,
00169 0.269266719309996,
00170 0.295524224714753,
00171 0.295524224714753,
00172 0.269266719309996,
00173 0.219086362515982,
00174 0.149451349150581,
00175 0.0666713443086881
00176 };
00177 static const G4double x[N] = {
00178 -0.973906528517172,
00179 -0.865063366688985,
00180 -0.679409568299024,
00181 -0.433395394129247,
00182 -0.148874338981631,
00183 0.148874338981631,
00184 0.433395394129247,
00185 0.679409568299024,
00186 0.865063366688985,
00187 0.973906528517172
00188 };
00189
00190 G4double Total = 0.0;
00191
00192
00193 for (G4int i = 0; i < N; i++)
00194 {
00195
00196 G4double KineticE = ((Up-Low)*x[i]+(Up+Low))/2.0;
00197
00198 Total += w[i]*ProbabilityDistributionFunction(fragment, KineticE);
00199
00200 }
00201 Total *= (Up-Low)/2.0;
00202 return Total;
00203 }
00204
00205
00207
00208
00209 G4double
00210 G4EvaporationProbability::ProbabilityDistributionFunction( const G4Fragment & fragment,
00211 G4double K)
00212 {
00213 G4int ResidualA = fragment.GetA_asInt() - theA;
00214 G4int ResidualZ = fragment.GetZ_asInt() - theZ;
00215 G4double U = fragment.GetExcitationEnergy();
00216
00217
00218
00219
00220
00221
00222
00223 G4double delta1 = fPairCorr->GetPairingCorrection(ResidualA,ResidualZ);
00224
00225 G4double delta0 = fPairCorr->GetPairingCorrection(fragment.GetA_asInt(),
00226 fragment.GetZ_asInt());
00227
00228
00229 G4double ParticleMass = fragment.ComputeGroundStateMass(theZ,theA);
00230 G4double ResidualMass = fragment.ComputeGroundStateMass(ResidualZ,ResidualA);
00231
00232 G4double theSeparationEnergy = ParticleMass + ResidualMass
00233 - fragment.GetGroundStateMass();
00234
00235 G4double a0 = theEvapLDPptr->LevelDensityParameter(fragment.GetA_asInt(),
00236 fragment.GetZ_asInt(),
00237 U - delta0);
00238
00239 G4double a1 = theEvapLDPptr->LevelDensityParameter(ResidualA, ResidualZ,
00240 U - theSeparationEnergy - delta1);
00241
00242
00243 G4double E0 = U - delta0;
00244
00245 G4double E1 = U - theSeparationEnergy - delta1 - K;
00246
00247 if (E1<0.) { return 0.; }
00248
00249
00250
00251
00252
00253
00254
00255 static const G4double pcoeff = millibarn/((pi*hbarc)*(pi*hbarc));
00256
00257
00258 G4double Prob = pcoeff*Gamma*ParticleMass*std::exp(2*(std::sqrt(a1*E1) - std::sqrt(a0*E0)))
00259 *K*CrossSection(fragment,K);
00260
00261 return Prob;
00262 }
00263
00264