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 #ifndef G4ProjectileFragmentCrossSection_h
00028 #define G4ProjectileFragmentCrossSection_h 1
00029
00030 #include <cmath>
00031 #include <iostream>
00032
00033
00034
00035
00036 class G4ProjectileFragmentCrossSection
00037 {
00038 public:
00039 G4ProjectileFragmentCrossSection()
00040 {
00041 p_S[1] = -2.38;
00042 p_S[2] = 0.27;
00043
00044 p_P[1] = -2.5840E+00;
00045 p_P[2] = -7.5700E-03;
00046
00047 p_Delta[1] = -1.0870E+00;
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;
00053 p_R[2] = -9.8160E-03;
00054
00055 p_Un[1] = 1.65;
00056
00057 p_Up[1] = 1.7880;
00058 p_Up[2] = +4.7210E-03;
00059 p_Up[3] = -1.3030E-05;
00060
00061 p_mn[1] = 0.400;
00062 p_mn[2] = 0.600;
00063
00064 p_mp[1] = -10.25;
00065 p_mp[2] = +10.1;
00066
00067 corr_d[1] = -25.0;
00068 corr_d[2] = 0.800;
00069 corr_r[1] = +20.0;
00070 corr_r[2] = 0.820;
00071 corr_y[1] = 200.0;
00072 corr_y[2] = 0.90;
00073 }
00074
00075 inline G4double doit(G4double Ap, G4double Zp, G4double At, G4double Zt, G4double A, G4double Z)
00076 {
00077
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
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
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
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
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
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
00131 zprob = zprob + 0.0020*A;
00132 cout <<"zprob "<<A<<" "<<dq<<" "<<Zp<<" "<<zbeta_p
00133 <<" "<<zbeta<<" "<<delta<<endl;
00134
00135
00136 G4double r = std::exp(p_R[1] + p_R[2]*A);
00137
00138
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
00147 if ((Zp-zbeta_p) < 0.0)
00148 {
00149 r=r*(1.0-0.0833*std::abs(Zp-zbeta_p));
00150 }
00151
00152
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
00157 G4double expo, fract;
00158 if((zprob-Z) > 0)
00159 {
00160
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
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
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 }
00185
00186 void testMe()
00187 {
00188 G4ProjectileFragmentCrossSection i;
00189 cout << i.doit(58, 28, 9, 4, 49, 28) << endl;
00190
00191 }
00192 private:
00193 G4double p_S[3];
00194 G4double p_P[3];
00195 G4double p_Delta[5];
00196 G4double p_R[3];
00197 G4double p_Un[2];
00198 G4double p_Up[4];
00199 G4double p_mn[3];
00200 G4double p_mp[3];
00201 G4double corr_d[3];
00202 G4double corr_r[3];
00203 G4double corr_y[3];
00204 };
00205 #endif