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