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 #include "G4ProtonEvaporationProbability.hh"
00038 #include "G4PhysicalConstants.hh"
00039 #include "G4SystemOfUnits.hh"
00040
00041 G4ProtonEvaporationProbability::G4ProtonEvaporationProbability() :
00042 G4EvaporationProbability(1,1,2,&theCoulombBarrier)
00043 {
00044 ResidualA = ResidualZ = theA = theZ = FragmentA = 0;
00045 ResidualAthrd = FragmentAthrd = U = 0.0;
00046 }
00047
00048 G4ProtonEvaporationProbability::~G4ProtonEvaporationProbability()
00049 {}
00050
00051 G4double G4ProtonEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)
00052 { return 1.0 + CCoeficient(fragment.GetZ_asInt()-GetZ());}
00053
00054 G4double G4ProtonEvaporationProbability::CalcBetaParam(const G4Fragment & )
00055 { return 0.0; }
00056
00057 G4double G4ProtonEvaporationProbability::CCoeficient(G4int aZ)
00058 {
00059
00060
00061
00062
00063
00064
00065
00066 G4double C = 0.0;
00067
00068 if (aZ >= 70) {
00069 C = 0.10;
00070 } else {
00071 C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
00072 }
00073
00074 return C;
00075
00076 }
00077
00079
00080
00081
00082
00083
00084
00085 G4double
00086 G4ProtonEvaporationProbability::CrossSection(const G4Fragment & fragment, G4double K)
00087 {
00088
00089
00090
00091 theA=GetA();
00092 theZ=GetZ();
00093 ResidualA=fragment.GetA_asInt()-theA;
00094 ResidualZ=fragment.GetZ_asInt()-theZ;
00095
00096 ResidualAthrd=fG4pow->Z13(ResidualA);
00097 FragmentA=fragment.GetA_asInt();
00098 FragmentAthrd=fG4pow->Z13(FragmentA);
00099
00100 U=fragment.GetExcitationEnergy();
00101
00102 if (OPTxs==0) {std::ostringstream errOs;
00103 errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (protons)!!" <<G4endl;
00104 throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00105 return 0.;}
00106 else if( OPTxs==1 ) return GetOpt1( K);
00107 else if( OPTxs==2 ||OPTxs==4) return GetOpt2( K);
00108 else if (OPTxs==3 ) return GetOpt3( K);
00109 else{
00110 std::ostringstream errOs;
00111 errOs << "BAD PROTON CROSS SECTION OPTION AT EVAPORATION!!" <<G4endl;
00112 throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00113 return 0.;
00114 }
00115 }
00116
00117
00118
00119
00120 G4double G4ProtonEvaporationProbability::GetOpt1(G4double K)
00121 {
00122 G4double Kc=K;
00123
00124
00125 if (K > 50*MeV) { Kc = 50*MeV; }
00126
00127 G4double landa, landa0, landa1, mu, mum0, mu1,nu, nu0, nu1, nu2,xs;
00128 G4double p, p0, p1, p2,Ec,delta,q,r,ji;
00129
00130 p0 = 15.72;
00131 p1 = 9.65;
00132 p2 = -449.0;
00133 landa0 = 0.00437;
00134 landa1 = -16.58;
00135 mum0 = 244.7;
00136 mu1 = 0.503;
00137 nu0 = 273.1;
00138 nu1 = -182.4;
00139 nu2 = -1.872;
00140 delta=0.;
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
00146 G4double resmu1 = fG4pow->powZ(ResidualA,mu1);
00147 mu = mum0*resmu1;
00148 nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
00149 q = landa - nu/(Ec*Ec) - 2*p*Ec;
00150 r = mu + 2*nu/Ec + p*(Ec*Ec);
00151
00152 ji=std::max(Kc,Ec);
00153 if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
00154 else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
00155 if (xs <0.0) {xs=0.0;}
00156
00157 return xs;
00158 }
00159
00160
00161
00162 G4double G4ProtonEvaporationProbability::GetOpt2(G4double K)
00163 {
00164
00165 G4double eekin,ekin,ff1,ff2,ff3,r0,fac,fac1,fac2,b0,xine_th(0);
00166
00167
00168
00169
00170
00171 if(!useSICB && K<=theCoulombBarrier.GetCoulombBarrier(ResidualA,ResidualZ,U))
00172 { return 0.0; }
00173
00174 eekin=K;
00175 G4int rnneu=ResidualA-ResidualZ;
00176 ekin=eekin/1000;
00177 r0=1.36*1.e-15;
00178 fac=pi*r0*r0;
00179 b0=2.247-0.915*(1.-1./ResidualAthrd);
00180 fac1=b0*(1.-1./ResidualAthrd);
00181 fac2=1.;
00182 if(rnneu > 1.5) { fac2 = fG4pow->logZ(rnneu); }
00183 xine_th= 1.e+31*fac*fac2*(1.+ResidualAthrd-fac1);
00184 xine_th=(1.-0.15*std::exp(-ekin))*xine_th/(1.00-0.0007*ResidualA);
00185 ff1=0.70-0.0020*ResidualA;
00186 ff2=1.00+1/G4double(ResidualA);
00187 ff3=0.8+18/G4double(ResidualA)-0.002*ResidualA;
00188 fac=1.-(1./(1.+std::exp(-8.*ff1*(std::log10(ekin)+1.37*ff2))));
00189 xine_th=xine_th*(1.+ff3*fac);
00190 ff1=1.-1/G4double(ResidualA)-0.001*ResidualA;
00191 ff2=1.17-2.7/G4double(ResidualA)-0.0014*ResidualA;
00192 fac=-8.*ff1*(std::log10(ekin)+2.0*ff2);
00193 fac=1./(1.+std::exp(fac));
00194 xine_th=xine_th*fac;
00195 if (xine_th < 0.0){
00196 std::ostringstream errOs;
00197 G4cout<<"WARNING: negative Wellisch cross section "<<G4endl;
00198 errOs << "RESIDUAL: A=" << ResidualA << " Z=" << ResidualZ <<G4endl;
00199 errOs <<" xsec("<<ekin<<" MeV) ="<<xine_th <<G4endl;
00200 throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00201 }
00202 return xine_th;
00203 }
00204
00205
00206 G4double G4ProtonEvaporationProbability::GetOpt3(const G4double K)
00207 {
00208
00209
00210
00211 G4double landa, landa0, landa1, mu, mum0, mu1,nu, nu0, nu1, nu2;
00212 G4double p, p0, p1, p2;
00213 p0 = 15.72;
00214 p1 = 9.65;
00215 p2 = -300.;
00216 landa0 = 0.00437;
00217 landa1 = -16.58;
00218 mum0 = 244.7;
00219 mu1 = 0.503;
00220 nu0 = 273.1;
00221 nu1 = -182.4;
00222 nu2 = -1.872;
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232 G4double ec,ecsq,xnulam,etest(0.),ra(0.),a,w,c,signor(1.),signor2,sig;
00233 G4double b,ecut,cut,ecut2,geom,elab;
00234
00235 G4double flow = 1.e-18;
00236 G4double spill= 1.e+18;
00237
00238 if (ResidualA <= 60) { signor = 0.92; }
00239 else if (ResidualA < 100) { signor = 0.8 + ResidualA*0.002; }
00240
00241 ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
00242 ecsq = ec * ec;
00243 p = p0 + p1/ec + p2/ecsq;
00244 landa = landa0*ResidualA + landa1;
00245 a = fG4pow->powZ(ResidualA,mu1);
00246 mu = mum0 * a;
00247 nu = a* (nu0+nu1*ec+nu2*ecsq);
00248
00249 c =std::min(3.15,ec*0.5);
00250 w = 0.7 * c / 3.15;
00251
00252 xnulam = nu / landa;
00253 if (xnulam > spill) { xnulam=0.; }
00254 if (xnulam >= flow) { etest =std::sqrt(xnulam) + 7.; }
00255
00256 a = -2.*p*ec + landa - nu/ecsq;
00257 b = p*ecsq + mu + 2.*nu/ec;
00258 ecut = 0.;
00259 cut = a*a - 4.*p*b;
00260 if (cut > 0.) { ecut = std::sqrt(cut); }
00261 ecut = (ecut-a) / (p+p);
00262 ecut2 = ecut;
00263
00264
00265
00266
00267 if (cut < 0.) { ecut2 = ecut; }
00268 elab = K * FragmentA /G4double(ResidualA);
00269 sig = 0.;
00270 if (elab <= ec) {
00271 if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
00272
00273 signor2 = (ec-elab-c) / w;
00274 signor2 = 1. + std::exp(signor2);
00275 sig = sig / signor2;
00276 }
00277 else{
00278 sig = (landa*elab+mu+nu/elab) * signor;
00279 geom = 0.;
00280
00281 if (xnulam < flow || elab < etest)
00282 {
00283 if (sig <0.0) {sig=0.0;}
00284 return sig;
00285 }
00286 geom = std::sqrt(theA*K);
00287 geom = 1.23*ResidualAthrd + ra + 4.573/geom;
00288 geom = 31.416 * geom * geom;
00289 sig = std::max(geom,sig);
00290
00291 }
00292 return sig;
00293 }
00294