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