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 "G4DeuteronEvaporationProbability.hh"
00039 #include "G4SystemOfUnits.hh"
00040
00041 G4DeuteronEvaporationProbability::G4DeuteronEvaporationProbability() :
00042 G4EvaporationProbability(2,1,3,&theCoulombBarrier)
00043 {
00044 ResidualA = ResidualZ = theA = theZ = FragmentA = 0;
00045 ResidualAthrd = FragmentAthrd = 0.0;
00046 }
00047
00048 G4DeuteronEvaporationProbability::~G4DeuteronEvaporationProbability()
00049 {}
00050
00051 G4double G4DeuteronEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)
00052 {
00053 return 1.0 + CCoeficient(fragment.GetZ_asInt()-GetZ());
00054 }
00055
00056 G4double G4DeuteronEvaporationProbability::CalcBetaParam(const G4Fragment & )
00057 {
00058 return 0.0;
00059 }
00060
00061 G4double G4DeuteronEvaporationProbability::CCoeficient(G4int aZ)
00062 {
00063
00064
00065
00066
00067
00068
00069
00070
00071 G4double C = 0.0;
00072
00073 if (aZ >= 70) {
00074 C = 0.10;
00075 } else {
00076 C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
00077 }
00078
00079 return C/2.0;
00080
00081 }
00082
00083
00085
00086
00087
00088
00089
00090 G4double
00091 G4DeuteronEvaporationProbability::CrossSection(const G4Fragment & fragment, G4double K)
00092 {
00093 theA=GetA();
00094 theZ=GetZ();
00095 ResidualA=fragment.GetA_asInt()-theA;
00096 ResidualZ=fragment.GetZ_asInt()-theZ;
00097
00098 ResidualAthrd=fG4pow->Z13(ResidualA);
00099 FragmentA=fragment.GetA_asInt();
00100 FragmentAthrd=fG4pow->Z13(FragmentA);
00101
00102 if (OPTxs==0) {std::ostringstream errOs;
00103 errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (deuterons)!!"
00104 <<G4endl;
00105 throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00106 return 0.;}
00107 if( OPTxs==1 || OPTxs==2) return G4DeuteronEvaporationProbability::GetOpt12( K);
00108 else if (OPTxs==3 || OPTxs==4) return G4DeuteronEvaporationProbability::GetOpt34( K);
00109 else{
00110 std::ostringstream errOs;
00111 errOs << "BAD Deuteron CROSS SECTION OPTION AT EVAPORATION!!" <<G4endl;
00112 throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00113 return 0.;
00114 }
00115 }
00116
00117
00118
00119
00120
00121 G4double G4DeuteronEvaporationProbability::GetOpt12(G4double K)
00122 {
00123 G4double Kc = K;
00124
00125
00126 if (K > 50*MeV) { Kc = 50*MeV; }
00127
00128 G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
00129
00130 G4double p0 = -38.21;
00131 G4double p1 = 922.6;
00132 G4double p2 = -2804.;
00133 G4double landa0 = -0.0323;
00134 G4double landa1 = -5.48;
00135 G4double mum0 = 336.1;
00136 G4double mu1 = 0.48;
00137 G4double nu0 = 524.3;
00138 G4double nu1 = -371.8;
00139 G4double nu2 = -5.924;
00140 G4double delta=1.2;
00141
00142 Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
00143 p = p0 + p1/Ec + p2/(Ec*Ec);
00144 landa = landa0*ResidualA + landa1;
00145 G4double resmu1 = fG4pow->powZ(ResidualA,mu1);
00146 mu = mum0*resmu1;
00147 nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
00148 q = landa - nu/(Ec*Ec) - 2*p*Ec;
00149 r = mu + 2*nu/Ec + p*(Ec*Ec);
00150
00151 ji=std::max(Kc,Ec);
00152 if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
00153 else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
00154
00155 if (xs <0.0) {xs=0.0;}
00156
00157 return xs;
00158 }
00159
00160
00161 G4double G4DeuteronEvaporationProbability::GetOpt34(G4double K)
00162
00163 {
00164
00165 G4double landa, mu, nu, p ,signor(1.),sig;
00166 G4double ec,ecsq,xnulam,etest(0.),a;
00167 G4double b,ecut,cut,ecut2,geom,elab;
00168
00169 G4double flow = 1.e-18;
00170 G4double spill= 1.e+18;
00171
00172 G4double p0 = 0.798;
00173 G4double p1 = 420.3;
00174 G4double p2 = -1651.;
00175 G4double landa0 = 0.00619;
00176 G4double landa1 = -7.54;
00177 G4double mum0 = 583.5;
00178 G4double mu1 = 0.337;
00179 G4double nu0 = 421.8;
00180 G4double nu1 = -474.5;
00181 G4double nu2 = -3.592;
00182
00183 G4double ra=0.80;
00184
00185
00186
00187 ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
00188 ecsq = ec * ec;
00189 p = p0 + p1/ec + p2/ecsq;
00190 landa = landa0*ResidualA + landa1;
00191 a = fG4pow->powZ(ResidualA,mu1);
00192 mu = mum0 * a;
00193 nu = a* (nu0+nu1*ec+nu2*ecsq);
00194 xnulam = nu / landa;
00195 if (xnulam > spill) { xnulam=0.; }
00196 if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); }
00197
00198 a = -2.*p*ec + landa - nu/ecsq;
00199 b = p*ecsq + mu + 2.*nu/ec;
00200 ecut = 0.;
00201 cut = a*a - 4.*p*b;
00202 if (cut > 0.) { ecut = std::sqrt(cut); }
00203 ecut = (ecut-a) / (p+p);
00204 ecut2 = ecut;
00205
00206
00207
00208
00209 if (cut < 0.) { ecut2 = ecut; }
00210 elab = K * FragmentA / G4double(ResidualA);
00211 sig = 0.;
00212
00213 if (elab <= ec) {
00214 if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
00215 }
00216 else {
00217 sig = (landa*elab+mu+nu/elab) * signor;
00218 geom = 0.;
00219 if (xnulam < flow || elab < etest) { return sig; }
00220 geom = std::sqrt(theA*K);
00221 geom = 1.23*ResidualAthrd + ra + 4.573/geom;
00222 geom = 31.416 * geom * geom;
00223 sig = std::max(geom,sig);
00224 }
00225 return sig;
00226 }
00227