27 #ifndef G4ProjectileFragmentCrossSection_h
28 #define G4ProjectileFragmentCrossSection_h 1
47 p_Delta[1] = -1.0870E+00;
48 p_Delta[2] = +3.0470E-02;
49 p_Delta[3] = +2.1353E-04;
50 p_Delta[4] = +7.1350E+01;
58 p_Up[2] = +4.7210E-03;
59 p_Up[3] = -1.3030E-05;
80 G4double S = p_S[2] * (At13 + Ap13 + p_S[1]);
82 G4double p = std::exp(p_P[2]*Ap + p_P[1]);
83 G4double yield_a = p * S * std::exp(-p * (Ap - A));
84 cout <<
"debug1 "<<yield_a<<endl;
89 f_mod_y=corr_y[1]*std::pow(A/Ap-corr_y[2], 2) + 1.0;
91 yield_a= yield_a * f_mod_y;
92 cout <<
"debug1 "<<yield_a<<endl;
95 G4double zbeta = A/(1.98+0.0155*std::pow(A, (2./3.)));
96 G4double zbeta_p = Ap/(1.98+0.0155*std::pow(Ap, (2./3.)));
100 delta = p_Delta[1] + p_Delta[2]*A;
104 delta = p_Delta[3]*A*A;
111 f_mod = corr_d[1]*std::pow(A/Ap-corr_d[2], 2) + 1.0;
121 cout <<
"dq "<<A<<
" "<<Ap<<
" "<<p_mp[1]
122 <<
" "<<p_mp[2]<<
" "<<dq<<
" "<<p_mp[1] + A/Ap*p_mp[2]<<endl;
126 dq = p_mn[1]*std::pow(A/Ap, 2.0) + p_mn[2]*std::pow(A/Ap, 4.0);
128 zprob = zprob + dq * (Zp-zbeta_p);
131 zprob = zprob + 0.0020*A;
132 cout <<
"zprob "<<A<<
" "<<dq<<
" "<<Zp<<
" "<<zbeta_p
133 <<
" "<<zbeta<<
" "<<delta<<endl;
136 G4double r = std::exp(p_R[1] + p_R[2]*A);
140 if (A/Ap > corr_r[2])
142 f_mod = corr_r[1]*Ap*std::pow(A/Ap-corr_r[2], 4.0)+1.0;
147 if ((Zp-zbeta_p) < 0.0)
149 r=r*(1.0-0.0833*std::abs(Zp-zbeta_p));
154 G4double u_p = p_Up[1] + p_Up[2]*A + p_Up[3]*A*A;
161 expo = -r*std::pow(std::abs(zprob-Z), u_n);
162 fract = std::exp(expo)*std::sqrt(r/3.14159);
167 expo = -r*std::pow(std::abs(zprob-Z), u_p);
168 fract = std::exp(expo)*std::sqrt(r/3.14159);
169 cout <<
"1 "<<expo<<
" "<<r<<
" "<<zprob<<
" "<<Z<<
" "<<u_p<<endl;
171 G4double dfdz = 1.2 + 0.647*std::pow(A/2.,0.3);
172 G4double z_exp = zprob + dfdz * std::log(10.) / (2.*
r);
175 expo = -r*std::pow(std::abs(zprob-z_exp), u_p);
176 fract = std::exp(expo)*std::sqrt(r/3.14159)
177 / std::pow(std::pow(10, dfdz), Z-z_exp);
181 cout <<
"debug "<<fract<<
" "<<yield_a<<endl;
189 cout << i.
doit(58, 28, 9, 4, 49, 28) << endl;
G4double doit(G4double Ap, G4double Zp, G4double At, G4double Zt, G4double A, G4double Z)
G4ProjectileFragmentCrossSection()