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 "G4PreCompoundDeuteron.hh"
00044 #include "G4SystemOfUnits.hh"
00045 #include "G4Deuteron.hh"
00046
00047 G4PreCompoundDeuteron::G4PreCompoundDeuteron()
00048 : G4PreCompoundIon(G4Deuteron::Deuteron(), &theDeuteronCoulombBarrier)
00049 {
00050 theA = GetA();
00051 theZ = GetZ();
00052 ResidualA = ResidualZ = 0;
00053 ResidualAthrd = FragmentAthrd = 0.0;
00054 FragmentA = theA;
00055 }
00056
00057 G4PreCompoundDeuteron::~G4PreCompoundDeuteron()
00058 {}
00059
00060 G4double G4PreCompoundDeuteron::FactorialFactor(G4int N, G4int P)
00061 {
00062 return G4double((N-1)*(N-2)*(P-1)*P)/2.0;
00063 }
00064
00065 G4double G4PreCompoundDeuteron::CoalescenceFactor(G4int A)
00066 {
00067 return 16.0/G4double(A);
00068 }
00069
00070 G4double G4PreCompoundDeuteron::GetRj(G4int nParticles, G4int nCharged)
00071 {
00072 G4double rj = 0.0;
00073 if(nCharged >=1 && (nParticles-nCharged) >=1) {
00074 G4double denominator = G4double(nParticles*(nParticles-1));
00075 rj = 2*nCharged*(nParticles-nCharged)/denominator;
00076 }
00077 return rj;
00078 }
00079
00081
00082
00083
00084
00085
00086 G4double G4PreCompoundDeuteron::CrossSection(G4double K)
00087 {
00088 ResidualA = GetRestA();
00089 ResidualZ = GetRestZ();
00090 theA = GetA();
00091 theZ = GetZ();
00092 ResidualAthrd = ResidualA13();
00093 FragmentA = theA + ResidualA;
00094 FragmentAthrd = g4pow->Z13(FragmentA);
00095
00096 if (OPTxs==0) { return GetOpt0( K); }
00097 else if( OPTxs==1 || OPTxs==2) { return GetOpt12( K); }
00098 else if (OPTxs==3 || OPTxs==4) { return GetOpt34( K); }
00099 else{
00100 std::ostringstream errOs;
00101 errOs << "BAD DEUTERON CROSS SECTION OPTION !!" <<G4endl;
00102 throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00103 return 0.;
00104 }
00105 }
00106
00107 G4double G4PreCompoundDeuteron::GetAlpha()
00108 {
00109 G4double C = 0.0;
00110 G4int aZ = theZ + ResidualZ;
00111 if (aZ >= 70)
00112 {
00113 C = 0.10;
00114 }
00115 else
00116 {
00117 C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
00118 }
00119 return 1.0 + C/2.0;
00120 }
00121
00122
00123
00124
00125 G4double G4PreCompoundDeuteron::GetOpt12(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 = -38.21;
00135 G4double p1 = 922.6;
00136 G4double p2 = -2804.;
00137 G4double landa0 = -0.0323;
00138 G4double landa1 = -5.48;
00139 G4double mm0 = 336.1;
00140 G4double mu1 = 0.48;
00141 G4double nu0 = 524.3;
00142 G4double nu1 = -371.8;
00143 G4double nu2 = -5.924;
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 G4double resmu1 = g4pow->powZ(ResidualA,mu1);
00150 mu = mm0*resmu1;
00151 nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
00152 q = landa - nu/(Ec*Ec) - 2*p*Ec;
00153 r = mu + 2*nu/Ec + p*(Ec*Ec);
00154
00155 ji=std::max(Kc,Ec);
00156 if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
00157 else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
00158
00159 if (xs <0.0) {xs=0.0;}
00160
00161 return xs;
00162 }
00163
00164
00165 G4double G4PreCompoundDeuteron::GetOpt34(G4double K)
00166
00167 {
00168
00169 G4double landa, mu, nu, p ,signor(1.),sig;
00170 G4double ec,ecsq,xnulam,etest(0.),a;
00171 G4double b,ecut,cut,ecut2,geom,elab;
00172
00173 G4double flow = 1.e-18;
00174 G4double spill= 1.e+18;
00175
00176 G4double p0 = 0.798;
00177 G4double p1 = 420.3;
00178 G4double p2 = -1651.;
00179 G4double landa0 = 0.00619;
00180 G4double landa1 = -7.54;
00181 G4double mm0 = 583.5;
00182 G4double mu1 = 0.337;
00183 G4double nu0 = 421.8;
00184 G4double nu1 = -474.5;
00185 G4double nu2 = -3.592;
00186
00187 G4double ra=0.80;
00188
00189
00190
00191 ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
00192 ecsq = ec * ec;
00193 p = p0 + p1/ec + p2/ecsq;
00194 landa = landa0*ResidualA + landa1;
00195 a = g4pow->powZ(ResidualA,mu1);
00196 mu = mm0 * a;
00197 nu = a* (nu0+nu1*ec+nu2*ecsq);
00198 xnulam = nu / landa;
00199 if (xnulam > spill) { xnulam=0.; }
00200 if (xnulam >= flow) { etest = 1.2 *std::sqrt(xnulam); }
00201
00202 a = -2.*p*ec + landa - nu/ecsq;
00203 b = p*ecsq + mu + 2.*nu/ec;
00204 ecut = 0.;
00205 cut = a*a - 4.*p*b;
00206 if (cut > 0.) { ecut = std::sqrt(cut); }
00207 ecut = (ecut-a) / (p+p);
00208 ecut2 = ecut;
00209
00210
00211
00212
00213 if (cut < 0.) { ecut2 = ecut; }
00214 elab = K * FragmentA / G4double(ResidualA);
00215 sig = 0.;
00216
00217 if (elab <= ec) {
00218 if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
00219 }
00220 else {
00221 sig = (landa*elab+mu+nu/elab) * signor;
00222 geom = 0.;
00223 if (xnulam < flow || elab < etest) { return sig; }
00224 geom = std::sqrt(theA*K);
00225 geom = 1.23*ResidualAthrd + ra + 4.573/geom;
00226 geom = 31.416 * geom * geom;
00227 sig = std::max(geom,sig);
00228 }
00229 return sig;
00230 }