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 "G4NeutronEvaporationProbability.hh"
00038 #include "G4SystemOfUnits.hh"
00039
00040 G4NeutronEvaporationProbability::G4NeutronEvaporationProbability() :
00041 G4EvaporationProbability(1,0,2,&theCoulombBarrier)
00042 {
00043 ResidualA = ResidualZ = theA = theZ = FragmentA = 0;
00044 ResidualAthrd = FragmentAthrd = 0.0;
00045 }
00046
00047 G4NeutronEvaporationProbability::~G4NeutronEvaporationProbability()
00048 {}
00049
00050 G4double G4NeutronEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment)
00051 {
00052 return 0.76+2.2/fG4pow->Z13(fragment.GetA_asInt()-GetA());
00053 }
00054
00055 G4double G4NeutronEvaporationProbability::CalcBetaParam(const G4Fragment & fragment)
00056 {
00057 return (2.12/fG4pow->Z23(fragment.GetA_asInt()-GetA()) - 0.05)*MeV/
00058 CalcAlphaParam(fragment);
00059 }
00060
00062
00063
00064
00065
00066
00067 G4double
00068 G4NeutronEvaporationProbability::CrossSection(const G4Fragment & fragment, G4double K)
00069 {
00070 theA=GetA();
00071 theZ=GetZ();
00072 ResidualA=fragment.GetA_asInt()-theA;
00073 ResidualZ=fragment.GetZ_asInt()-theZ;
00074
00075 ResidualAthrd=fG4pow->Z13(ResidualA);
00076 FragmentA=fragment.GetA_asInt();
00077 FragmentAthrd=fG4pow->Z13(FragmentA);
00078
00079 if (OPTxs==0) {std::ostringstream errOs;
00080 errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (neutrons)!!" <<G4endl;
00081 throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00082 return 0.;}
00083 else if( OPTxs==1 ||OPTxs==2) return GetOpt12( K);
00084 else if (OPTxs==3 || OPTxs==4) return GetOpt34( K);
00085 else{
00086 std::ostringstream errOs;
00087 errOs << "BAD NEUTRON CROSS SECTION OPTION AT EVAPORATION!!" <<G4endl;
00088 throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00089 return 0.;
00090 }
00091 }
00092
00093
00094
00095
00096 G4double G4NeutronEvaporationProbability::GetOpt12(G4double K)
00097 {
00098 G4double Kc=K;
00099
00100
00101
00102
00103 if (K > 50*MeV) { Kc = 50*MeV; }
00104
00105 G4double landa, landa0, landa1, mu, mum0, mu1,nu, nu0, nu1, nu2,xs;
00106
00107 landa0 = 18.57;
00108 landa1 = -22.93;
00109 mum0 = 381.7;
00110 mu1 = 24.31;
00111 nu0 = 0.172;
00112 nu1 = -15.39;
00113 nu2 = 804.8;
00114 landa = landa0/ResidualAthrd + landa1;
00115 mu = mum0*ResidualAthrd + mu1*ResidualAthrd*ResidualAthrd;
00116 nu = nu0*ResidualAthrd*ResidualA + nu1*ResidualAthrd*ResidualAthrd + nu2 ;
00117 xs=landa*Kc + mu + nu/Kc;
00118 if (xs <= 0.0 ){
00119 std::ostringstream errOs;
00120 G4cout<<"WARNING: NEGATIVE OPT=1 neutron cross section "<<G4endl;
00121 errOs << "RESIDUAL: Ar=" << ResidualA << " Zr=" << ResidualZ <<G4endl;
00122 errOs <<" xsec("<<Kc<<" MeV) ="<<xs <<G4endl;
00123 throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00124 }
00125 return xs;
00126 }
00127
00128
00129 G4double G4NeutronEvaporationProbability::GetOpt34(G4double K)
00130 {
00131 G4double landa, landa0, landa1, mu, mum0, mu1,nu, nu0, nu1, nu2;
00132 G4double p, p0;
00133 G4double flow,ec,ecsq,xnulam,etest(0.),ra(0.),a,signor(1.),sig;
00134 G4double b,ecut,cut,ecut2,geom,elab;
00135
00136
00137 landa0=0;
00138 landa1=0;
00139 mum0=0.;
00140 mu1=0.;
00141 nu0=0.;
00142 nu1=0.;
00143 nu2=0.;
00144 p0=0.;
00145
00146 flow = 1.e-18;
00147
00148
00149 p0 = -312.;
00150 landa0 = 12.10;
00151 landa1= -11.27;
00152 mum0 = 234.1;
00153 mu1 = 38.26;
00154 nu0 = 1.55;
00155 nu1 = -106.1;
00156 nu2 = 1280.8;
00157
00158 if (ResidualA < 40) { signor =0.7 + ResidualA*0.0075; }
00159 if (ResidualA > 210) { signor = 1. + (ResidualA-210)/250.; }
00160 landa = landa0/ResidualAthrd + landa1;
00161 mu = mum0*ResidualAthrd + mu1*ResidualAthrd*ResidualAthrd;
00162 nu = nu0*ResidualAthrd*ResidualA + nu1*ResidualAthrd*ResidualAthrd + nu2;
00163
00164
00165 if (nu < 0.) { nu=-nu; }
00166
00167 ec = 0.5;
00168 ecsq = 0.25;
00169 p = p0;
00170 xnulam = 1.;
00171 etest = 32.;
00172
00173
00174
00175
00176 a = -2.*p*ec + landa - nu/ecsq;
00177 b = p*ecsq + mu + 2.*nu/ec;
00178 ecut = 0.;
00179 cut = a*a - 4.*p*b;
00180 if (cut > 0.) { ecut = std::sqrt(cut); }
00181 ecut = (ecut-a) / (p+p);
00182 ecut2 = ecut;
00183 if (cut < 0.) { ecut2 = ecut - 2.; }
00184 elab = K * FragmentA / G4double(ResidualA);
00185 sig = 0.;
00186 if (elab <= ec) {
00187 if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
00188 }
00189 else {
00190 sig = (landa*elab+mu+nu/elab) * signor;
00191 geom = 0.;
00192 if (xnulam < flow || elab < etest) { return sig; }
00193 geom = std::sqrt(theA*K);
00194 geom = 1.23*ResidualAthrd + ra + 4.573/geom;
00195 geom = 31.416 * geom * geom;
00196 sig = std::max(geom,sig);
00197
00198 }
00199 return sig;
00200 }
00201