76 {
77
80 G4double S = p_S[2] * (At13 + Ap13 + p_S[1]);
81
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;
85
87 if (A/Ap > corr_y[2])
88 {
89 f_mod_y=corr_y[1]*std::pow(A/Ap-corr_y[2], 2) + 1.0;
90 }
91 yield_a= yield_a * f_mod_y;
92 cout << "debug1 "<<yield_a<<endl;
93
94
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.)));
98 if(A > p_Delta[4])
99 {
100 delta = p_Delta[1] + p_Delta[2]*A;
101 }
102 else
103 {
104 delta = p_Delta[3]*A*A;
105 }
106
107
109 if(A/Ap > corr_d[2])
110 {
111 f_mod = corr_d[1]*std::pow(A/Ap-corr_d[2], 2) + 1.0;
112 }
113 delta = delta*f_mod;
115
116
118 if((Zp-zbeta_p)>0)
119 {
121 cout << "dq "<<A<<" "<<Ap<<" "<<p_mp[1]
122 <<" "<<p_mp[2]<<" "<<dq<<" "<<p_mp[1] + A/Ap*p_mp[2]<<endl;
123 }
124 else
125 {
126 dq = p_mn[1]*std::pow(A/Ap, 2.0) + p_mn[2]*std::pow(A/Ap, 4.0);
127 }
128 zprob = zprob + dq * (Zp-zbeta_p);
129
130
131 zprob = zprob + 0.0020*A;
132 cout <<"zprob "<<A<<" "<<dq<<" "<<Zp<<" "<<zbeta_p
133 <<" "<<zbeta<<" "<<delta<<endl;
134
135
136 G4double r = std::exp(p_R[1] + p_R[2]*A);
137
138
139 f_mod=1.0;
140 if (A/Ap > corr_r[2])
141 {
142 f_mod = corr_r[1]*Ap*std::pow(A/Ap-corr_r[2], 4.0)+1.0;
143 }
144 r = r*f_mod;
145
146
147 if ((Zp-zbeta_p) < 0.0)
148 {
149 r=r*(1.0-0.0833*std::abs(Zp-zbeta_p));
150 }
151
152
154 G4double u_p = p_Up[1] + p_Up[2]*A + p_Up[3]*A*A;
155
156
158 if((zprob-Z) > 0)
159 {
160
161 expo = -r*std::pow(std::abs(zprob-Z), u_n);
162 fract = std::exp(expo)*std::sqrt(r/3.14159);
163 }
164 else
165 {
166
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;
170
171 G4double dfdz = 1.2 + 0.647*std::pow(A/2.,0.3);
172 G4double z_exp = zprob + dfdz * std::log(10.) / (2.*r);
173 if( Z>z_exp )
174 {
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);
178 }
179 }
180
181 cout << "debug "<<fract<<" "<<yield_a<<endl;
183 return epaxv2;
184 }