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 "G4PreCompoundTriton.hh"
00044 #include "G4SystemOfUnits.hh"
00045 #include "G4Triton.hh"
00046
00047 G4PreCompoundTriton::G4PreCompoundTriton()
00048 : G4PreCompoundIon(G4Triton::Triton(), &theTritonCoulombBarrier)
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 G4PreCompoundTriton::~G4PreCompoundTriton()
00060 {}
00061
00062 G4double G4PreCompoundTriton::FactorialFactor(G4int N, const G4int P)
00063 {
00064 return G4double((N-3)*(P-2)*(N-2)*(P-1)*(N-1)*P)/6.0;
00065 }
00066
00067 G4double G4PreCompoundTriton::CoalescenceFactor(G4int A)
00068 {
00069 return 243.0/G4double(A*A);
00070 }
00071
00072 G4double G4PreCompoundTriton::GetRj(G4int nParticles, G4int nCharged)
00073 {
00074 G4double rj = 0.0;
00075 if(nCharged >= 1 && (nParticles-nCharged) >= 2) {
00076 G4double denominator =
00077 G4double(nParticles*(nParticles-1)*(nParticles-2));
00078 rj = G4double(3*nCharged*(nParticles-nCharged)*(nParticles-nCharged-1))
00079 /denominator;
00080 }
00081 return rj;
00082 }
00083
00085
00086
00087
00088
00089
00090 G4double G4PreCompoundTriton::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 TRITON CROSS SECTION OPTION !!" <<G4endl;
00106 throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00107 return 0.;
00108 }
00109 }
00110
00111 G4double G4PreCompoundTriton::GetAlpha()
00112 {
00113 G4double C = 0.0;
00114 G4int aZ = theZ + ResidualZ;
00115 if (aZ >= 70)
00116 {
00117 C = 0.10;
00118 }
00119 else
00120 {
00121 C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
00122 }
00123
00124 return 1.0 + C/3.0;
00125 }
00126
00127
00128
00129
00130
00131 G4double G4PreCompoundTriton::GetOpt12(G4double K)
00132 {
00133 G4double Kc=K;
00134
00135
00136 if (K > 50*MeV) { Kc=50*MeV; }
00137
00138 G4double landa ,mu ,nu ,p , Ec,q,r,ji,xs;
00139
00140 G4double p0 = -11.04;
00141 G4double p1 = 619.1;
00142 G4double p2 = -2147.;
00143 G4double landa0 = -0.0426;
00144 G4double landa1 = -10.33;
00145 G4double mm0 = 601.9;
00146 G4double mu1 = 0.37;
00147 G4double nu0 = 583.0;
00148 G4double nu1 = -546.2;
00149 G4double nu2 = 1.718;
00150 G4double delta=1.2;
00151
00152 Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
00153 p = p0 + p1/Ec + p2/(Ec*Ec);
00154 landa = landa0*ResidualA + landa1;
00155
00156 G4double resmu1 = g4pow->powZ(ResidualA,mu1);
00157 mu = mm0*resmu1;
00158 nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
00159 q = landa - nu/(Ec*Ec) - 2*p*Ec;
00160 r = mu + 2*nu/Ec + p*(Ec*Ec);
00161
00162 ji=std::max(Kc,Ec);
00163 if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
00164 else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
00165
00166 if (xs <0.0) {xs=0.0;}
00167
00168 return xs;
00169 }
00170
00171
00172 G4double G4PreCompoundTriton::GetOpt34(G4double K)
00173
00174 {
00175 G4double landa, mu, nu, p , signor(1.),sig;
00176 G4double ec,ecsq,xnulam,etest(0.),a;
00177 G4double b,ecut,cut,ecut2,geom,elab;
00178
00179 G4double flow = 1.e-18;
00180 G4double spill= 1.e+18;
00181
00182 G4double p0 = -21.45;
00183 G4double p1 = 484.7;
00184 G4double p2 = -1608.;
00185 G4double landa0 = 0.0186;
00186 G4double landa1 = -8.90;
00187 G4double mm0 = 686.3;
00188 G4double mu1 = 0.325;
00189 G4double nu0 = 368.9;
00190 G4double nu1 = -522.2;
00191 G4double nu2 = -4.998;
00192
00193 G4double ra=0.80;
00194
00195
00196
00197 ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
00198 ecsq = ec * ec;
00199 p = p0 + p1/ec + p2/ecsq;
00200 landa = landa0*ResidualA + landa1;
00201 a = g4pow->powZ(ResidualA,mu1);
00202 mu = mm0 * a;
00203 nu = a* (nu0+nu1*ec+nu2*ecsq);
00204 xnulam = nu / landa;
00205 if (xnulam > spill) { xnulam=0.; }
00206 if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); }
00207
00208 a = -2.*p*ec + landa - nu/ecsq;
00209 b = p*ecsq + mu + 2.*nu/ec;
00210 ecut = 0.;
00211 cut = a*a - 4.*p*b;
00212 if (cut > 0.) { ecut = std::sqrt(cut); }
00213 ecut = (ecut-a) / (p+p);
00214 ecut2 = ecut;
00215
00216
00217
00218
00219 if (cut < 0.) { ecut2 = ecut; }
00220 elab = K * FragmentA / G4double(ResidualA);
00221 sig = 0.;
00222
00223 if (elab <= ec) {
00224 if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
00225 }
00226 else {
00227 sig = (landa*elab+mu+nu/elab) * signor;
00228 geom = 0.;
00229 if (xnulam < flow || elab < etest) { return sig; }
00230 geom = std::sqrt(theA*K);
00231 geom = 1.23*ResidualAthrd + ra + 4.573/geom;
00232 geom = 31.416 * geom * geom;
00233 sig = std::max(geom,sig);
00234 }
00235 return sig;
00236 }