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