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