#include <G4ProjectileFragmentCrossSection.hh>
Public Member Functions | |
G4ProjectileFragmentCrossSection () | |
G4double | doit (G4double Ap, G4double Zp, G4double At, G4double Zt, G4double A, G4double Z) |
void | testMe () |
Definition at line 36 of file G4ProjectileFragmentCrossSection.hh.
G4ProjectileFragmentCrossSection::G4ProjectileFragmentCrossSection | ( | ) | [inline] |
Definition at line 39 of file G4ProjectileFragmentCrossSection.hh.
00040 { 00041 p_S[1] = -2.38; // scale factor for xsect in barn 00042 p_S[2] = 0.27; 00043 00044 p_P[1] = -2.5840E+00; // slope of mass yield curve 00045 p_P[2] = -7.5700E-03; 00046 00047 p_Delta[1] = -1.0870E+00; // centroid rel. to beta-stability 00048 p_Delta[2] = +3.0470E-02; 00049 p_Delta[3] = +2.1353E-04; 00050 p_Delta[4] = +7.1350E+01; 00051 00052 p_R[1] = +0.885E+00; // width parameter R 00053 p_R[2] = -9.8160E-03; 00054 00055 p_Un[1] = 1.65; // slope par. n-rich ride of Z distr. 00056 00057 p_Up[1] = 1.7880; // slope par. p-rich ride of Z distr. 00058 p_Up[2] = +4.7210E-03; 00059 p_Up[3] = -1.3030E-05; 00060 00061 p_mn[1] = 0.400; // memory effect n-rich projectiles 00062 p_mn[2] = 0.600; 00063 00064 p_mp[1] = -10.25; // memory effect p-rich projectiles 00065 p_mp[2] = +10.1; 00066 00067 corr_d[1] = -25.0; // correction close to proj.: centroid dzp 00068 corr_d[2] = 0.800; 00069 corr_r[1] = +20.0; // correction close to proj.: width R 00070 corr_r[2] = 0.820; 00071 corr_y[1] = 200.0; // correction close to proj.: Yield_a 00072 corr_y[2] = 0.90; 00073 }
G4double G4ProjectileFragmentCrossSection::doit | ( | G4double | Ap, | |
G4double | Zp, | |||
G4double | At, | |||
G4double | Zt, | |||
G4double | A, | |||
G4double | Z | |||
) | [inline] |
Definition at line 75 of file G4ProjectileFragmentCrossSection.hh.
Referenced by testMe().
00076 { 00077 // calculate mass yield 00078 G4double Ap13 = std::pow(Ap, 1./3.); 00079 G4double At13 = std::pow(At, 1./3.); 00080 G4double S = p_S[2] * (At13 + Ap13 + p_S[1]); 00081 // cout << "debug0 "<<S<<" "<<At13<<" "<<Ap13<<" "<<p_S[1]<<" "<<p_S[2]<<endl; 00082 G4double p = std::exp(p_P[2]*Ap + p_P[1]); 00083 G4double yield_a = p * S * std::exp(-p * (Ap - A)); 00084 cout << "debug1 "<<yield_a<<endl; 00085 // modification close to projectile 00086 G4double f_mod_y=1.0; 00087 if (A/Ap > corr_y[2]) 00088 { 00089 f_mod_y=corr_y[1]*std::pow(A/Ap-corr_y[2], 2) + 1.0; 00090 } 00091 yield_a= yield_a * f_mod_y; 00092 cout << "debug1 "<<yield_a<<endl; 00093 00094 // calculate maximum of charge dispersion zprob 00095 G4double zbeta = A/(1.98+0.0155*std::pow(A, (2./3.))); 00096 G4double zbeta_p = Ap/(1.98+0.0155*std::pow(Ap, (2./3.))); 00097 G4double delta; 00098 if(A > p_Delta[4]) 00099 { 00100 delta = p_Delta[1] + p_Delta[2]*A; 00101 } 00102 else 00103 { 00104 delta = p_Delta[3]*A*A; 00105 } 00106 00107 // modification close to projectile 00108 G4double f_mod=1.0; 00109 if(A/Ap > corr_d[2]) 00110 { 00111 f_mod = corr_d[1]*std::pow(A/Ap-corr_d[2], 2) + 1.0; 00112 } 00113 delta = delta*f_mod; 00114 G4double zprob = zbeta+delta; 00115 00116 // correction for proton- and neutron-rich projectiles 00117 G4double dq; 00118 if((Zp-zbeta_p)>0) 00119 { 00120 dq = std::exp(p_mp[1] + G4double(A)/G4double(Ap)*p_mp[2]); 00121 cout << "dq "<<A<<" "<<Ap<<" "<<p_mp[1] 00122 <<" "<<p_mp[2]<<" "<<dq<<" "<<p_mp[1] + A/Ap*p_mp[2]<<endl; 00123 } 00124 else 00125 { 00126 dq = p_mn[1]*std::pow(A/Ap, 2.0) + p_mn[2]*std::pow(A/Ap, 4.0); 00127 } 00128 zprob = zprob + dq * (Zp-zbeta_p); 00129 00130 // small corr. since Xe-129 and Pb-208 are not on Z_beta line 00131 zprob = zprob + 0.0020*A; 00132 cout <<"zprob "<<A<<" "<<dq<<" "<<Zp<<" "<<zbeta_p 00133 <<" "<<zbeta<<" "<<delta<<endl; 00134 00135 // calculate width parameter R 00136 G4double r = std::exp(p_R[1] + p_R[2]*A); 00137 00138 // modification close to projectile 00139 f_mod=1.0; 00140 if (A/Ap > corr_r[2]) 00141 { 00142 f_mod = corr_r[1]*Ap*std::pow(A/Ap-corr_r[2], 4.0)+1.0; 00143 } 00144 r = r*f_mod; 00145 00146 // change width according to dev. from beta-stability 00147 if ((Zp-zbeta_p) < 0.0) 00148 { 00149 r=r*(1.0-0.0833*std::abs(Zp-zbeta_p)); 00150 } 00151 00152 // calculate slope parameters u_n, u_p 00153 G4double u_n = p_Un[1]; 00154 G4double u_p = p_Up[1] + p_Up[2]*A + p_Up[3]*A*A; 00155 00156 // calculate charge dispersion 00157 G4double expo, fract; 00158 if((zprob-Z) > 0) 00159 { 00160 // neutron-rich 00161 expo = -r*std::pow(std::abs(zprob-Z), u_n); 00162 fract = std::exp(expo)*std::sqrt(r/3.14159); 00163 } 00164 else 00165 { 00166 // proton-rich 00167 expo = -r*std::pow(std::abs(zprob-Z), u_p); 00168 fract = std::exp(expo)*std::sqrt(r/3.14159); 00169 cout << "1 "<<expo<<" "<<r<<" "<<zprob<<" "<<Z<<" "<<u_p<<endl; 00170 // go to exponential slope 00171 G4double dfdz = 1.2 + 0.647*std::pow(A/2.,0.3); 00172 G4double z_exp = zprob + dfdz * std::log(10.) / (2.*r); 00173 if( Z>z_exp ) 00174 { 00175 expo = -r*std::pow(std::abs(zprob-z_exp), u_p); 00176 fract = std::exp(expo)*std::sqrt(r/3.14159) 00177 / std::pow(std::pow(10, dfdz), Z-z_exp); 00178 } 00179 } 00180 00181 cout << "debug "<<fract<<" "<<yield_a<<endl; 00182 G4double epaxv2=fract*yield_a; 00183 return epaxv2; 00184 }
void G4ProjectileFragmentCrossSection::testMe | ( | ) | [inline] |
Definition at line 186 of file G4ProjectileFragmentCrossSection.hh.
References doit().
00187 { 00188 G4ProjectileFragmentCrossSection i; 00189 cout << i.doit(58, 28, 9, 4, 49, 28) << endl; 00190 // Sigma = 9.800163E-13 b 00191 }