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