Heavy ion Analysis Libriares
Loading...
Searching...
No Matches
FemtoFsiParsed.cxx
1/*
2 * FemtoFsiParsed.cxx
3 *
4 * Created on: 30 kwi 2018
5 * Author: Daniel Wielanek
6 * E-mail: daniel.wielanek@gmail.com
7 * Warsaw University of Technology, Faculty of Physics
8 */
9#include "FemtoFsiParsed.h"
10
11#include <TMath.h>
12#include <iostream>
13
14namespace Hal {
15 FemtoFsiParsed::FemtoFsiParsed() {
16 fFsiAapi = {.26, .25, .348, -.0637, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
17 -.028, -.082, .2795, -.0086, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
18 fFsiApin = {1.2265, 15.63, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
19 -.75, 76.88, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
20 fFsiAand = {-3.295, -881.9, 4537., 28645., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
21 -13.837, 11.505, 0., 10.416, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
22 -32.18, 10.14, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
23 -60.213, 13.33, 0., -70.309, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
24 fFsiAadd = {1061.7, .003194, 568.49, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
25 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
26 -1085., 19870., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
27 fFsiAakk = {.247677, .01947977, .2997015, .9604, .96511, .094, .11, .333, .222};
28 fsi_aapap__ = {-.94, -1.98, .1, -1.4, .37, .1, .3, .267, -.01, 1.66, .553, -.01};
29 fsi_fdh__ = {
30 17., 7.8, 23.7, 2230.1218, .225, .081, -.063, -.65, -2.73, .137, -.071, -.148,
31 .112, 1e-6, 1e-6, -.36, 1e-6, 1e-6, 1.344, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6,
32 1e-6, -5.628, 2.18, 2.4, 2.81, .001, -5.4, -5.4, -5.4, 0., 0., 0.,
33 0., -6.35, -11.88, 0., 0., 0., 0., 0., 0., 0., 0., 0.,
34 0., 0., 0., 0., 0., 0., 0., 0., 1.93, 1.84, .5, .001,
35 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
36 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
37 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
38 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
39 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
40 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
41 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
42 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
43 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
44 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
45 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
46 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
47 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
48 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
49 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
50 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
51 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
52 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
53 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
54 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
55 2.7, 2.8, 2.7, 1.12139906, -44.36, 64., 784.9, 477.9, 2.27, 0., 0., 0.,
56 0., 0., 0., 0., 0., 0., -69.973, 0., 0., 0., 0., 0.,
57 0., 3.529, 3.19, 3.15, 2.95, 0., 1.7, 1.7, 1.7, 0., 0., 0.,
58 0., 2., 2.63, 0., 0., 0., 0., 0., 0., 0., 0., 0.,
59 0., 0., 0., 0., 0., 0., 0., 0., 3.35, 3.37, 2.95, 0.,
60 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
61 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
62 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
63 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
64 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
65 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
66 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
67 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
68 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
69 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
70 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
71 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
72 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
73 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
74 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
75 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
76 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
77 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
78 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
79 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
80 .1149517, .1046257, .1148757, .118601, .7947389, 2.281208, 8.7, .4, .1561219, 0., 0., 0.,
81 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
82 0., 0., .1013293, .1020966, .1080476, .001, .1847221, .1847221, .1847221, 0., 0., 0.,
83 0., .4, .1150687, 0., 0., 0., 0., 0., 0., 0., 0., 0.,
84 0., 0., 0., 0., 0., 0., 0., 0., .09736083, .0970831, .1080476, .001,
85 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
86 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
87 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
88 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
89 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
90 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
91 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
92 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
93 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
94 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
95 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
96 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
97 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
98 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
99 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
100 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
101 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
102 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
103 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
104 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
105 2.545739, 2.779789, 2.585795, 5.023544, .124673, .392518, .09, 2., 4.058058, 0., 0., 0.,
106 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
107 0., 0., 2.252623, 2.278575, 2.234089, .001, 2.003144, 2.003144, 2.003144, 0., 0., 0.,
108 0., 2., 4.132163, 0., 0., 0., 0., 0., 0., 0., 0., 0.,
109 0., 0., 0., 0., 0., 0., 0., 0., 2.272703, 2.256355, 2.234089, .001,
110 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
111 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
112 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
113 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
114 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
115 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
116 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
117 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
118 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
119 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
120 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
121 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
122 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
123 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
124 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
125 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
126 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
127 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
128 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
129 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
130 fsi_rhoh__ = {.25, .25, .25, 1., 1., 1., 1., .3333, .3333, 1., 1., 1., 1., 1., 1., 1., 1., .1111, 1., .25,
131 1., 1., 1., .3333, .25, 1., .25, .25, .25, .25, .75, .75, .75, 0., 0., 0., 0., .6667, .6667, 0.,
132 0., 0., 0., 0., 0., 0., 0., .3333, 0., .75, 0., 0., 0., .6667, .75, 0., .75, .75, .75, .75,
133 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .5556, 0., 0.,
134 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
135 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
136 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
137 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
138 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
139 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
140 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
141 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
142 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
143 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
144 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
145 fsi_amch__ = {
146 .93956563, .93827231, .93956563, 3.72737978, .13957, .13498, .13957, .93956563, .93827231, .13957,
147 .13957, .13957, .13957, .493677, .493677, .493677, .493677, 1.87561339, 1.87561339, 2.80892165,
148 2.80892165, .497672, .497672, 1.87561339, .93827231, .93827231, .93827231, .93956563, 1.115684, .93827231,
149 .93956563, .93827231, .93827231, 3.72737978, .13957, .13498, .13957, 1.87561339, 1.87561339, .493677,
150 .493677, .93827231, .93827231, .493677, .493677, .93827231, .93827231, 1.87561339, 3.72737978, 2.80892165,
151 3.72737978, .497672, .497672, 2.80892165, 2.80892165, 3.72737978, 1.115684, 1.115684, 1.115684, .93827231,
152 0., 1., 0., 2., 1., 0., 1., 0., 1., 1.,
153 1., 1., -1., 1., 1., 1., -1., 1., 1., 1.,
154 1., 0., 0., 1., 1., 1., 1., 0., 0., 1.,
155 0., 1., 1., 2., -1., 0., 1., 1., 1., -1.,
156 1., 1., 1., -1., 1., 1., 1., 1., 2., 1.,
157 2., 0., 0., 1., 1., 2., 0., 0., 0., -1.,
158 2, 2, 2, 1, 1, 1, 1, 2, 2, 1,
159 1, 1, 1, 1, 1, 1, 1, 3, 1, 2,
160 1, 1, 1, 2, 2, 1, 2, 2, 2, 2};
161 c__1 = 1;
162 c_b5 = {(Double_t) 1., (Double_t) 0.};
163 c_b2 = {1., 0.};
164 c__2 = 2;
165 c_b33 = {.6667, 0.};
166 c_b34 = {.3333, 0.};
167 ifirst = 0;
168 fEta = 0;
169 fIch = 0;
170 fJrat = 0;
171 }
172
173 void FemtoFsiParsed::Fsiw(int j, Double_t& weif, Double_t& wei, Double_t& wein) {
174 Double_t c1n, c2n;
175 c1n = fFsiPoc.c1 * fFsiPoc.cn;
176 if (c1n != 0.) { fFsiPoc.ac1 = 137.036 / (c1n * fFsiMom.e1); }
177 /* m1-->E1 */
178 c2n = fFsiPoc.c2 * fFsiPoc.cn;
179 if (c2n != 0.) { fFsiPoc.ac2 = 137.036 / (c2n * fFsiMom.e2); }
180 /* ----------------------------------------------------------- */
181 /* m2-->E2 */
182
183 Fsipn(weif);
184
185 /* weight due to particle-nucleus FSI */
186 fJrat = 0;
187 Boosttoprf();
188 Fsiwf(wei);
189
190 /* weight due to particle-particle-nucleus FSI */
191 wein = wei;
192 if (fFsiNs.i3c * j != 0) {
193 fFsiFppn.ff12.r = 1., fFsiFppn.ff12.i = 0.;
194 fFsiFppn.ff21.r = 1., fFsiFppn.ff21.i = 0.;
195 fJrat = 1;
196 Vz(wein);
197
198 /* weight due to particle-particle FSI */
199 }
200 }
201
202 void FemtoFsiParsed::Ltran(Double_t* p0, Double_t* p, Double_t* ps) {
203 Double_t d1, d2, d3;
204
205 /* Local variables */
206 static Double_t h, am0, pp0, p0s, epm;
207
208 /* ==>calculating particle 4-momentum PS={PSX,PSY,PSZ,ES} */
209 /* in rest frame of a system 0 with 4-momentum P0={P0X,P0Y,P0Z,E0} */
210 /* from its 4-momentum P={PX,PY,PZ,E} */
211 /* ----------------------------------------------------------------------- */
212 /* Parameter adjustments */
213 --ps;
214 --p;
215 --p0;
216
217 /* Function Body */
218 /* Computing 2nd power */
219 d1 = p0[1];
220 /* Computing 2nd power */
221 d2 = p0[2];
222 /* Computing 2nd power */
223 d3 = p0[3];
224 p0s = d1 * d1 + d2 * d2 + d3 * d3;
225 /* Computing 2nd power */
226 d1 = p0[4];
227 am0 = TMath::Sqrt(d1 * d1 - p0s);
228 epm = p0[4] + am0;
229 pp0 = p[1] * p0[1] + p[2] * p0[2] + p[3] * p0[3];
230 h = (pp0 / epm - p[4]) / am0;
231 ps[1] = p[1] + p0[1] * h;
232 ps[2] = p[2] + p0[2] * h;
233 ps[3] = p[3] + p0[3] * h;
234 ps[4] = (p0[4] * p[4] - pp0) / am0;
235 }
236
237 void FemtoFsiParsed::Ltranb(Double_t* p0, Double_t* ps, Double_t* p) {
238 Double_t d__1, d__2, d__3;
239
240 /* Local variables */
241 Double_t hs, am0, p0s, epm, psp0;
242
243 /* ==>calculating particle 4-momentum P={PX,PY,PZ,E} */
244 /* from its 4-momentum PS={PSX,PSY,PSZ,ES} */
245 /* in rest frame of a system 0 with 4-momentum P0={P0X,P0Y,P0Z,E0} */
246 /* ----------------------------------------------------------------------- */
247 /* Parameter adjustments */
248 --p;
249 --ps;
250 --p0;
251
252 /* Function Body */
253 /* Computing 2nd power */
254 d__1 = p0[1];
255 /* Computing 2nd power */
256 d__2 = p0[2];
257 /* Computing 2nd power */
258 d__3 = p0[3];
259 p0s = d__1 * d__1 + d__2 * d__2 + d__3 * d__3;
260 /* Computing 2nd power */
261 d__1 = p0[4];
262 am0 = TMath::Sqrt(d__1 * d__1 - p0s);
263 epm = p0[4] + am0;
264 psp0 = ps[1] * p0[1] + ps[2] * p0[2] + ps[3] * p0[3];
265 hs = (psp0 / epm + ps[4]) / am0;
266 p[1] = ps[1] + p0[1] * hs;
267 p[2] = ps[2] + p0[2] * hs;
268 p[3] = ps[3] + p0[3] * hs;
269 p[4] = (p0[4] * ps[4] + psp0) / am0;
270 }
271
272 void FemtoFsiParsed::Ltran12() {
273 Double_t d__1, d__2, d__3;
274 /* Local variables */
275 Double_t h1, ee, p112, p1s, p2s, p12s;
276
277 /* ==>calculating particle momentum in PRF {EE,PPX,PPY,PPZ} from */
278 /* - the momentum of the first particle {E1,P1X,P1Y,P1Z) in NRF */
279 /* part. momenta in NRF */
280 /* 4-coord. of emis. points in NRF */
281 /* fm --> 1/GeV */
282 fFsiCoor.x1 *= fFsiCons.w;
283 fFsiCoor.y1 *= fFsiCons.w;
284 fFsiCoor.z1 *= fFsiCons.w;
285 fFsiCoor.t1 *= fFsiCons.w;
286 fFsiCoor.x2 *= fFsiCons.w;
287 fFsiCoor.y2 *= fFsiCons.w;
288 fFsiCoor.z2 *= fFsiCons.w;
289 fFsiCoor.t2 *= fFsiCons.w;
290 /* calculating Ri, Pi and Ei */
291 fFsiCoor.r1 = TMath::Sqrt(fFsiCoor.x1 * fFsiCoor.x1 + fFsiCoor.y1 * fFsiCoor.y1 + fFsiCoor.z1 * fFsiCoor.z1);
292 fFsiCoor.r2 = TMath::Sqrt(fFsiCoor.x2 * fFsiCoor.x2 + fFsiCoor.y2 * fFsiCoor.y2 + fFsiCoor.z2 * fFsiCoor.z2);
293 p1s = fFsiMom.p1x * fFsiMom.p1x + fFsiMom.p1y * fFsiMom.p1y + fFsiMom.p1z * fFsiMom.p1z;
294 p2s = fFsiMom.p2x * fFsiMom.p2x + fFsiMom.p2y * fFsiMom.p2y + fFsiMom.p2z * fFsiMom.p2z;
295 fFsiMom.p1 = TMath::Sqrt(p1s);
296 fFsiMom.p2 = TMath::Sqrt(p2s);
297 fFsiMom.e1 = TMath::Sqrt(fFsiPoc.am1 * fFsiPoc.am1 + p1s);
298 fFsiMom.e2 = TMath::Sqrt(fFsiPoc.am2 * fFsiPoc.am2 + p2s);
299 /* ----------------------------------------------------------------------- */
300 fSi_p12.e12 = fFsiMom.e1 + fFsiMom.e2;
301 fSi_p12.p12x = fFsiMom.p1x + fFsiMom.p2x;
302 fSi_p12.p12y = fFsiMom.p1y + fFsiMom.p2y;
303 fSi_p12.p12z = fFsiMom.p1z + fFsiMom.p2z;
304 /* Computing 2nd power */
305 d__1 = fSi_p12.p12x;
306 /* Computing 2nd power */
307 d__2 = fSi_p12.p12y;
308 /* Computing 2nd power */
309 d__3 = fSi_p12.p12z;
310 p12s = d__1 * d__1 + d__2 * d__2 + d__3 * d__3;
311 /* Computing 2nd power */
312 d__1 = fSi_p12.e12;
313 fSi_p12.am12 = TMath::Sqrt(d__1 * d__1 - p12s);
314 fSi_p12.epm = fSi_p12.e12 + fSi_p12.am12;
315 fSi_p12.p12 = TMath::Sqrt(p12s);
316 p112 = fFsiMom.p1x * fSi_p12.p12x + fFsiMom.p1y * fSi_p12.p12y + fFsiMom.p1z * fSi_p12.p12z;
317 h1 = (p112 / fSi_p12.epm - fFsiMom.e1) / fSi_p12.am12;
318 fFsiPrf.ppx = fFsiMom.p1x + fSi_p12.p12x * h1;
319 fFsiPrf.ppy = fFsiMom.p1y + fSi_p12.p12y * h1;
320 fFsiPrf.ppz = fFsiMom.p1z + fSi_p12.p12z * h1;
321 ee = (fSi_p12.e12 * fFsiMom.e1 - p112) / fSi_p12.am12;
322 /* Computing 2nd power */
323 d__1 = ee;
324 /* Computing 2nd power */
325 d__2 = fFsiPoc.am1;
326 fFsiPrf.aks = d__1 * d__1 - d__2 * d__2;
327 fFsiPrf.ak = TMath::Sqrt(fFsiPrf.aks);
328 /* W WRITE(6,38)'AK ',AK,'K ',PPX,PPY,PPZ,EE */
329 /* L38: */
330 }
331
332 void FemtoFsiParsed::Fsipn(Double_t& weif) {
333 /* System generated locals */
334 Double_t d__1, d__2;
335 complex z__1, z__2;
336 /* Local variables */
337 Double_t hs, xh;
338 Double_t rhos;
339
340 /* calculating particle-nucleus Coulomb Wave functions FFij */
341 /* part. momenta in NRF */
342 /* 4-coord. of emis. points in NRF */
343 fFsiFppn.ff12.r = 1., fFsiFppn.ff12.i = 0.;
344 fFsiFppn.ff21.r = 1., fFsiFppn.ff21.i = 0.;
345 /* ACH=1.D0 */
346 /* WEIF=1.D0 */
347 if (fFsiNs.i3c == 0) { return; }
348 fIch = (int) fFsiPoc.c1;
349 if (fIch == 0) { goto L11; }
350 xh = fFsiPoc.ac1 * fFsiMom.p1;
351 fFsiAch._1.ach = Acp(xh);
352 fFsiAch._1.achr = TMath::Sqrt(fFsiAch._1.ach);
353 fEta = 0.;
354 if (xh != 0.) { fEta = 1 / xh; }
355 rhos = fFsiMom.p1 * fFsiCoor.r1;
356 hs = fFsiCoor.x1 * fFsiMom.p1x + fFsiCoor.y1 * fFsiMom.p1y + fFsiCoor.z1 * fFsiMom.p1z;
357 Ff1(z__2, rhos, hs);
358 z__1.r = fFsiFppn.ff12.r * z__2.r - fFsiFppn.ff12.i * z__2.i, z__1.i = fFsiFppn.ff12.r * z__2.i + fFsiFppn.ff12.i * z__2.r;
359 fFsiFppn.ff12.r = z__1.r, fFsiFppn.ff12.i = z__1.i;
360 if (fFsiNs.iqs == 0) { goto L11; }
361 rhos = fFsiMom.p1 * fFsiCoor.r2;
362 hs = fFsiCoor.x2 * fFsiMom.p1x + fFsiCoor.y2 * fFsiMom.p1y + fFsiCoor.z2 * fFsiMom.p1z;
363 Ff1(z__2, rhos, hs);
364 z__1.r = fFsiFppn.ff21.r * z__2.r - fFsiFppn.ff21.i * z__2.i, z__1.i = fFsiFppn.ff21.r * z__2.i + fFsiFppn.ff21.i * z__2.r;
365 fFsiFppn.ff21.r = z__1.r, fFsiFppn.ff21.i = z__1.i;
366 L11:
367 fIch = (int) fFsiPoc.c2;
368 if (fIch == 0) { goto L10; }
369 xh = fFsiPoc.ac2 * fFsiMom.p2;
370 fFsiAch._1.ach = Acp(xh);
371 fFsiAch._1.achr = TMath::Sqrt(fFsiAch._1.ach);
372 fEta = 0.;
373 if (xh != 0.) { fEta = 1 / xh; }
374 rhos = fFsiMom.p2 * fFsiCoor.r2;
375 hs = fFsiCoor.x2 * fFsiMom.p2x + fFsiCoor.y2 * fFsiMom.p2y + fFsiCoor.z2 * fFsiMom.p2z;
376 Ff1(z__2, rhos, hs);
377 z__1.r = fFsiFppn.ff12.r * z__2.r - fFsiFppn.ff12.i * z__2.i, z__1.i = fFsiFppn.ff12.r * z__2.i + fFsiFppn.ff12.i * z__2.r;
378 fFsiFppn.ff12.r = z__1.r, fFsiFppn.ff12.i = z__1.i;
379 /* W WRITE(6,41)'AC2 ',AC2,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS */
380 /* L41: */
381 /* W WRITE(6,40)'FF12 ',DREAL(FF12),DIMAG(FF12) */
382 if (fFsiNs.iqs == 0) { goto L10; }
383 rhos = fFsiMom.p2 * fFsiCoor.r1;
384 hs = fFsiCoor.x1 * fFsiMom.p2x + fFsiCoor.y1 * fFsiMom.p2y + fFsiCoor.z1 * fFsiMom.p2z;
385 Ff1(z__2, rhos, hs);
386 z__1.r = fFsiFppn.ff21.r * z__2.r - fFsiFppn.ff21.i * z__2.i, z__1.i = fFsiFppn.ff21.r * z__2.i + fFsiFppn.ff21.i * z__2.r;
387 fFsiFppn.ff21.r = z__1.r, fFsiFppn.ff21.i = z__1.i;
388 /* W WRITE(6,41)'AC1 ',AC1,'ACH ',ACH,'ETA ',ETA,'RHOS ',RHOS,'HS ',HS */
389 /* W WRITE(6,40)'FF21 ',DREAL(FF21),DIMAG(FF21) */
390 /* L40: */
391 L10:
392 /* WEIF = the weight due to the Coulomb particle-nucleus interaction */
393 /* Computing 2nd power */
394 d__1 = fFsiFppn.ff12.r;
395 /* Computing 2nd power */
396 d__2 = d_imag(&fFsiFppn.ff12);
397 weif = d__1 * d__1 + d__2 * d__2;
398 if (fFsiNs.iqs == 1) {
399 /* Computing 2nd power */
400 d__1 = fFsiFppn.ff21.r;
401 /* Computing 2nd power */
402 d__2 = d_imag(&fFsiFppn.ff21);
403 weif = (weif + d__1 * d__1 + d__2 * d__2) * .5;
404 }
405 }
406
407 FemtoFsiParsed::~FemtoFsiParsed() {
408 // TODO Auto-generated destructor stub
409 }
410
411 Double_t FemtoFsiParsed::Gpipi(Double_t x, int j) {
412 /* System generated locals */
413 Double_t ret_val, d__1;
414
415 /* Local variables */
416 Double_t om, xx;
417
418 /* --- GPIPI = k*COTG(DELTA), X=k^2 */
419 /* -- J=1(2) corresponds to isospin=0(2) */
420 om = TMath::Sqrt(x + fsi_c__._1.ams);
421 xx = x / fsi_c__._1.ams;
422 ret_val = om / fFsiAapi.e_1[j * 20 - 20];
423 /* Computing 2nd power */
424 d__1 = fFsiAapi.e_1[j * 20 - 20];
425 ret_val *= (fFsiAapi.e_1[j * 20 - 18] - d__1 * d__1) * xx + 1 + fFsiAapi.e_1[j * 20 - 17] * xx * xx;
426 ret_val /= (fFsiAapi.e_1[j * 20 - 18] + fFsiAapi.e_1[j * 20 - 19] / fFsiAapi.e_1[j * 20 - 20]) * xx + 1;
427 return ret_val;
428 } /* gpipi_ */
429
430 Double_t FemtoFsiParsed::Gpin(Double_t x, int j) {
431 /* System generated locals */
432 Double_t ret_val;
433
434 /* --- GPIN = k*COTG(DELTA), X=k^2 */
435 /* -- J=1(2) corresponds to piN isospin=1/2(3/2) */
436 ret_val = 1 / fFsiApin.e_1[j * 20 - 20] + fFsiApin.e_1[j * 20 - 19] * .5 * x;
437 return ret_val;
438 } /* gpin_ */
439
440 Double_t FemtoFsiParsed::Gnd(Double_t x, int j) {
441 /* System generated locals */
442 Double_t ret_val;
443
444 /* Local variables */
445 int i__;
446 Double_t xx;
447
448 /* --- GND = k*COTG(DELTA), X=k^2 */
449 /* --- J=1(2) corresp. to nd(pd), S=1/2, */
450 /* --- J=3(4) corresp. to nd(pd), S=3/2 */
451 xx = x;
452 ret_val = 1 / fFsiAand.e_1[j * 20 - 20] + fFsiAand.e_1[j * 20 - 19] * .5 * x;
453 for (i__ = 4; i__ <= 4; ++i__) {
454 xx *= x;
455 /* L1: */
456 ret_val += fFsiAand.e_1[i__ + j * 20 - 21] * xx;
457 }
458 ret_val /= fFsiAand.e_1[j * 20 - 18] * x + 1;
459 return ret_val;
460 } /* gnd_ */
461
462 Double_t FemtoFsiParsed::Gdd(Double_t x, int j) {
463 /* System generated locals */
464 Double_t ret_val;
465 /* Local variables */
466 Double_t e, er, tand;
467
468 /* --- GDD = k*COTG(DELTA), X=k^2 */
469 /* --- J=1,2,3 corresp. to S=0,1,2 */
470 e = x / 2 / fsi_c__._2.am;
471 er = TMath::Sqrt(e);
472 if (j == 1) {
473 ret_val = er * (fFsiAadd.e_1[0] * exp(-e / fFsiAadd.e_1[1]) - fFsiAadd.e_1[2]);
474 ret_val /= fFsiCons.dr;
475 /* from degree to radian */
476 tand = TMath::Tan(ret_val);
477 if (tand == 0.) { tand = 1e-10; }
478 ret_val = TMath::Sqrt(x) / tand;
479 }
480 if (j == 2) { ret_val = 1e10; }
481 if (j == 3) {
482 ret_val = er * (fFsiAadd.e_1[40] + fFsiAadd.e_1[41] * e);
483 ret_val /= fFsiCons.dr;
484 /* from degree to radian */
485 tand = TMath::Tan(ret_val);
486 if (tand == 0.) { tand = 1e-10; }
487 ret_val = TMath::Sqrt(x) / tand;
488 }
489 return ret_val;
490 } /* gdd_ */
491
492 Double_t FemtoFsiParsed::DfrSin(Double_t x) {
493 /* Initialized data */
494
495 Double_t c0 = 1.2533141373155003;
496 int nb = 15;
497 int nc1 = 25;
498 int nc2 = 28;
499 Double_t b[16] = {.630414043145705392,
500 -.423445114057053335,
501 .376171726433436566,
502 -.162494891545095674,
503 .038222557786330087,
504 -.005645634771321909,
505 5.74549519768974e-4,
506 -4.287071532102e-5,
507 2.451207499233e-6,
508 -1.10988418409e-7,
509 4.082497317e-9,
510 -1.24498302e-10,
511 3.200484e-12,
512 -7.0324e-14,
513 1.336e-15,
514 -2.2e-17};
515 Double_t c1[26] = {.990560479373497549,
516 -.012183509831478997,
517 -.00248274282311306,
518 2.66094952647247e-4,
519 -1.079068987406e-6,
520 -4.883681753933e-6,
521 9.99055266368e-7,
522 -7.5092717372e-8,
523 -1.9079487573e-8,
524 9.090797293e-9,
525 -1.966236033e-9,
526 1.64772911e-10,
527 6.3079714e-11,
528 -3.6432219e-11,
529 1.053693e-11,
530 -1.716438e-12,
531 -1.07124e-13,
532 2.04099e-13,
533 -9.0064e-14,
534 2.5506e-14,
535 -4.036e-15,
536 -5.7e-16,
537 7.62e-16,
538 -3.63e-16,
539 1.18e-16,
540 -2.5e-17};
541 Double_t c2[29] = {.0465577987375164561,
542 .044992130201239414,
543 -.0017542871396514532,
544 -1.465340025810678e-4,
545 3.91330408630159e-5,
546 -3.4932286597731e-6,
547 -3.153530032345e-7,
548 1.876582008529e-7,
549 -3.7755280493e-8,
550 2.665516501e-9,
551 1.0881448122e-9,
552 -5.355007671e-10,
553 1.315765447e-10,
554 -1.52860881e-11,
555 -3.3947646e-12,
556 2.7020267e-12,
557 -9.463142e-13,
558 2.071565e-13,
559 -1.26931e-14,
560 -1.39756e-14,
561 8.5929e-15,
562 -3.107e-15,
563 7.515e-16,
564 -6.48e-17,
565 -5.22e-17,
566 3.86e-17,
567 -1.65e-17,
568 5e-18,
569 -9e-19};
570
571 /* System generated locals */
572 Double_t ret_val, d__1;
573
574 /* Local variables */
575 Double_t h__;
576 int i__;
577 Double_t r__, s, v, y, b0, b1, b2, alfa;
578
579 v = TMath::Abs(x);
580 if (v < 8.) {
581 y = v * .125;
582 /* Computing 2nd power */
583 d__1 = y;
584 h__ = d__1 * d__1 * 2 - 1;
585 alfa = h__ + h__;
586 b1 = 0.;
587 b2 = 0.;
588 for (i__ = nb; i__ >= 0; --i__) {
589 b0 = b[i__] + alfa * b1 - b2;
590 b2 = b1;
591 /* L4: */
592 b1 = b0;
593 }
594 h__ = TMath::Sqrt(v) * y * (b0 - b2);
595 } else {
596 r__ = 1 / v;
597 h__ = r__ * 10 - 1;
598 alfa = h__ + h__;
599 b1 = 0.;
600 b2 = 0.;
601 for (i__ = nc1; i__ >= 0; --i__) {
602 b0 = c1[i__] + alfa * b1 - b2;
603 b2 = b1;
604 /* L5: */
605 b1 = b0;
606 }
607 s = b0 - h__ * b2;
608 b1 = 0.;
609 b2 = 0.;
610 for (i__ = nc2; i__ >= 0; --i__) {
611 b0 = c2[i__] + alfa * b1 - b2;
612 b2 = b1;
613 /* L6: */
614 b1 = b0;
615 }
616 h__ = c0 - TMath::Sqrt(r__) * (s * TMath::Cos(v) + (b0 - h__ * b2) * TMath::Sin(v));
617 }
618 if (x < 0.) { h__ = -h__; }
619 ret_val = h__;
620 return ret_val;
621 } /* dfrsin_ */
622
623 Double_t FemtoFsiParsed::DfrCos(Double_t x) {
624 /* Initialized data */
625
626 Double_t c0 = 1.2533141373155003;
627 int na = 16;
628 int nc1 = 25;
629 int nc2 = 28;
630 Double_t a[17] = {.764351386641860002,
631 -.431355475476601793,
632 .432881999797266531,
633 -.26973310338387111,
634 .084160453208769354,
635 -.01546524484461382,
636 .00187855423439822,
637 -1.62649776188875e-4,
638 1.0573976563833e-5,
639 -5.36093398892e-7,
640 2.1816584549e-8,
641 -7.29016212e-10,
642 2.0373325e-11,
643 -4.8344e-13,
644 9.865e-15,
645 -1.75e-16,
646 3e-18};
647 Double_t c1[26] = {.990560479373497549,
648 -.012183509831478997,
649 -.00248274282311306,
650 2.66094952647247e-4,
651 -1.079068987406e-6,
652 -4.883681753933e-6,
653 9.99055266368e-7,
654 -7.5092717372e-8,
655 -1.9079487573e-8,
656 9.090797293e-9,
657 -1.966236033e-9,
658 1.64772911e-10,
659 6.3079714e-11,
660 -3.6432219e-11,
661 1.053693e-11,
662 -1.716438e-12,
663 -1.07124e-13,
664 2.04099e-13,
665 -9.0064e-14,
666 2.5506e-14,
667 -4.036e-15,
668 -5.7e-16,
669 7.62e-16,
670 -3.63e-16,
671 1.18e-16,
672 -2.5e-17};
673 Double_t c2[29] = {.0465577987375164561,
674 .044992130201239414,
675 -.0017542871396514532,
676 -1.465340025810678e-4,
677 3.91330408630159e-5,
678 -3.4932286597731e-6,
679 -3.153530032345e-7,
680 1.876582008529e-7,
681 -3.7755280493e-8,
682 2.665516501e-9,
683 1.0881448122e-9,
684 -5.355007671e-10,
685 1.315765447e-10,
686 -1.52860881e-11,
687 -3.3947646e-12,
688 2.7020267e-12,
689 -9.463142e-13,
690 2.071565e-13,
691 -1.26931e-14,
692 -1.39756e-14,
693 8.5929e-15,
694 -3.107e-15,
695 7.515e-16,
696 -6.48e-17,
697 -5.22e-17,
698 3.86e-17,
699 -1.65e-17,
700 5e-18,
701 -9e-19};
702
703 /* System generated locals */
704 Double_t ret_val, d__1;
705
706 /* Local variables */
707 Double_t h__;
708 int i__;
709 Double_t r__, s, v, b0, b1, b2, alfa;
710
711 v = TMath::Abs(x);
712 if (v < 8.) {
713 /* Computing 2nd power */
714 d__1 = v;
715 h__ = d__1 * d__1 * .03125 - 1;
716 alfa = h__ + h__;
717 b1 = 0.;
718 b2 = 0.;
719 for (i__ = na; i__ >= 0; --i__) {
720 b0 = a[i__] + alfa * b1 - b2;
721 b2 = b1;
722 /* L1: */
723 b1 = b0;
724 }
725 h__ = TMath::Sqrt(v) * (b0 - h__ * b2);
726 } else {
727 r__ = 1 / v;
728 h__ = r__ * 10 - 1;
729 alfa = h__ + h__;
730 b1 = 0.;
731 b2 = 0.;
732 for (i__ = nc1; i__ >= 0; --i__) {
733 b0 = c1[i__] + alfa * b1 - b2;
734 b2 = b1;
735 /* L2: */
736 b1 = b0;
737 }
738 s = b0 - h__ * b2;
739 b1 = 0.;
740 b2 = 0.;
741 for (i__ = nc2; i__ >= 0; --i__) {
742 b0 = c2[i__] + alfa * b1 - b2;
743 b2 = b1;
744 /* L3: */
745 b1 = b0;
746 }
747 h__ = c0 - TMath::Sqrt(r__) * ((b0 - h__ * b2) * TMath::Cos(v) - s * TMath::Sin(v));
748 }
749 if (x < 0.) { h__ = -h__; }
750 ret_val = h__;
751 return ret_val;
752 } /* dfrcos_ */
753
754 void FemtoFsiParsed::Cgamma(complex& ret_val, complex& z__) {
755 /* Initialized data */
756
757 Double_t pi = (Double_t) 3.14159265358979324;
758 Double_t c1 = (Double_t) 2.5066282746310005;
759 Double_t c__[16] = {(Double_t) 41.624436916439068,
760 (Double_t) -51.224241022374774,
761 (Double_t) 11.338755813488977,
762 (Double_t) -.747732687772388,
763 (Double_t) .008782877493061,
764 (Double_t) -1.899030264e-6,
765 (Double_t) 1.946335e-9,
766 (Double_t) -1.99345e-10,
767 (Double_t) 8.433e-12,
768 (Double_t) 1.486e-12,
769 (Double_t) -8.06e-13,
770 (Double_t) 2.93e-13,
771 (Double_t) -1.02e-13,
772 (Double_t) 3.7e-14,
773 (Double_t) -1.4e-14,
774 (Double_t) 6e-15};
775
776 /* Format strings */
777 // static char fmt_101[] = "(\002ARGUMENT EQUALS NON-POSITIVE INTEGER =
778 // \\002,1p,e15.1)";
779 /* System generated locals */
780 int i__1;
781 complex q__1, q__2, q__3, q__4, q__5, q__6, q__7;
782
783 /* Local variables */
784 complex f, h__;
785 int k;
786 complex s, u, v;
787 Double_t x;
788 // char errtxt[80];
789
790 /* Fortran I/O blocks */
791 // icilist io___9 = { 0, errtxt, 0, fmt_101, 80, 1 };
792 u.r = z__.r, u.i = z__.i;
793 x = u.r;
794 if (r_imag(&u) == (Double_t) 0. && -TMath::Abs(x) == (Double_t)((int) x)) {
795 f.r = (Double_t) 0., f.i = (Double_t) 0.;
796 h__.r = (Double_t) 0., h__.i = (Double_t) 0.;
797
798 /* CALL MTLPRT(NAME,'C305.1',ERRTXT) */
799 } else {
800 if (x >= (Double_t) 1.) {
801 f.r = (Double_t) 1., f.i = (Double_t) 0.;
802 v.r = u.r, v.i = u.i;
803 } else if (x >= (Double_t) 0.) {
804 c_div(&q__1, &c_b5, &u);
805 f.r = q__1.r, f.i = q__1.i;
806 q__1.r = u.r + 1, q__1.i = u.i;
807 v.r = q__1.r, v.i = q__1.i;
808 } else {
809 f.r = (Double_t) 1., f.i = (Double_t) 0.;
810 q__1.r = 1 - u.r, q__1.i = -u.i;
811 v.r = q__1.r, v.i = q__1.i;
812 }
813 h__.r = (Double_t) 1., h__.i = (Double_t) 0.;
814 s.r = c__[0], s.i = (Double_t) 0.;
815 for (k = 1; k <= 15; ++k) {
816 q__3.r = v.r - k, q__3.i = v.i;
817 i__1 = k - 1;
818 q__4.r = v.r + i__1, q__4.i = v.i;
819 c_div(&q__2, &q__3, &q__4);
820 q__1.r = q__2.r * h__.r - q__2.i * h__.i, q__1.i = q__2.r * h__.i + q__2.i * h__.r;
821 h__.r = q__1.r, h__.i = q__1.i;
822 // std::cout<<"L"<<k<<" ("<<h__.r<<" "<<h__.i<<")"<<std::endl;
823 /* L1: */
824 i__1 = k;
825 q__2.r = c__[i__1] * h__.r, q__2.i = c__[i__1] * h__.i;
826 q__1.r = s.r + q__2.r, q__1.i = s.i + q__2.i;
827 s.r = q__1.r, s.i = q__1.i;
828 }
829 q__1.r = v.r + (Double_t) 4.5, q__1.i = v.i;
830 h__.r = q__1.r, h__.i = q__1.i;
831 q__6.r = v.r - (Double_t) .5, q__6.i = v.i;
832 c_log(&q__7, &h__);
833 q__5.r = q__6.r * q__7.r - q__6.i * q__7.i, q__5.i = q__6.r * q__7.i + q__6.i * q__7.r;
834 q__4.r = q__5.r - h__.r, q__4.i = q__5.i - h__.i;
835 c_exp(&q__3, &q__4);
836 q__2.r = c1 * q__3.r, q__2.i = c1 * q__3.i;
837 q__1.r = q__2.r * s.r - q__2.i * s.i, q__1.i = q__2.r * s.i + q__2.i * s.r;
838 h__.r = q__1.r, h__.i = q__1.i;
839 if (x < (Double_t) 0.) {
840 q__2.r = pi, q__2.i = (Double_t) 0.;
841 q__5.r = pi * u.r, q__5.i = pi * u.i;
842 c_sin(&q__4, &q__5);
843 q__3.r = q__4.r * h__.r - q__4.i * h__.i, q__3.i = q__4.r * h__.i + q__4.i * h__.r;
844 c_div(&q__1, &q__2, &q__3);
845 h__.r = q__1.r, h__.i = q__1.i;
846 }
847 }
848
849 q__1.r = f.r * h__.r - f.i * h__.i, q__1.i = f.r * h__.i + f.i * h__.r;
850 ret_val.r = q__1.r, ret_val.i = q__1.i;
851 z__.r = ret_val.r;
852 z__.i = ret_val.i;
853 return;
854 }
855
856 void FemtoFsiParsed::Ckkb() {
857 /*
858 Double_t S4 = fFsiPrf.aks + fFsiAakk.e_1[0];
859 Double_t S = 4*S4;
860 Double_t AKPIPI = TMath::Sqrt(S4 - fFsiAakk.e_1[1]);
861 Double_t AAKK_2 = fFsiAakk.e_1[2];
862 Double_t EETA2 = (S+fFsiAakk.e_1[2]-fFsiAakk.e_1[1])*
863 (S+fFsiAakk.e_1[2]-fFsiAakk.e_1[1])/4/S;
864 Double_t AKPIETA = TMath::Sqrt(EETA2 - fFsiAakk.e_1[2]);
865 SetC(TComplex(), 0);
866 */
867
868 /* System generated locals */
869 Double_t d__1, d__2, d__3;
870 complex z__1, z__2, z__3, z__4;
871
872 /* Local variables */
873 Double_t s, s4, eeta2, akpipi, akpieta;
874
875 /* saturated by S*(980) and delta(982) resonances */
876 /* calculates KK-b scattering amplitude, */
877 s4 = fFsiPrf.aks + fFsiAakk.e_1[0];
878 s = s4 * 4;
879 akpipi = TMath::Sqrt(s4 - fFsiAakk.e_1[1]);
880 /* Computing 2nd power */
881 d__1 = s + fFsiAakk.e_1[2] - fFsiAakk.e_1[1];
882 eeta2 = d__1 * d__1 / 4 / s;
883 akpieta = TMath::Sqrt(eeta2 - fFsiAakk.e_1[2]);
884 d__1 = fFsiAakk.e_1[5] / 2;
885 z__2.r = d__1, z__2.i = 0.;
886 d__2 = fFsiAakk.e_1[3] - s;
887 d__3 = -fFsiPrf.ak * fFsiAakk.e_1[5] - akpipi * fFsiAakk.e_1[6];
888 z__3.r = d__2, z__3.i = d__3;
889 z_div(&z__1, &z__2, &z__3);
890 fsi_c__._2.c__[0].r = z__1.r, fsi_c__._2.c__[0].i = z__1.i;
891 d__1 = fFsiAakk.e_1[7] / 2;
892 z__3.r = d__1, z__3.i = 0.;
893 d__2 = fFsiAakk.e_1[4] - s;
894 d__3 = -fFsiPrf.ak * fFsiAakk.e_1[7] - akpieta * fFsiAakk.e_1[8];
895 z__4.r = d__2, z__4.i = d__3;
896 z_div(&z__2, &z__3, &z__4);
897 z__1.r = fsi_c__._2.c__[0].r + z__2.r, z__1.i = fsi_c__._2.c__[0].i + z__2.i;
898 fsi_c__._2.c__[0].r = z__1.r, fsi_c__._2.c__[0].i = z__1.i;
899
900 } /* ckkb_ */
901
902 void FemtoFsiParsed::Cpap() {
903 Double_t AM2 = .93956563;
904 ;
905 Double_t AM = fsi_c__._2.am;
906 Double_t AMU2_AMU1 = AM2 / AM;
907 Double_t AKS = fFsiPrf.aks;
908 Double_t AK2S0 = AMU2_AMU1 * AKS;
909 Double_t AK2S = fsi_2cha__.ak2s;
910 AK2S = AK2S0 - 2 * AM2 * (AM2 - AM);
911 // TComplex AK2(fsi_2cha__.ak2.r,fsi_2cha__.ak2.i);
912 if (AK2S >= 0) {
913 // AK2= TComplex(TMath::Sqrt(AK2S),0);//k2
914 fsi_2cha__.ak2 = TMath::Sqrt(AK2S);
915 } else {
916 fsi_2cha__.ak2 = TMath::Sqrt(-AK2S);
917 // AK2 = TComplex(0,TMath::Sqrt(-AK2S)); //kappa2
918 }
919 Int_t ISPIN = fFsiAch._1.ispin;
920 Double_t a_re = fsi_aapap__1.aapapr[fFsiAch._1.ispin * 3 - 1] * fFsiPrf.aks / 2 - .016 - fsi_2cha__.hcp2;
921 Double_t a_im = fsi_aapap__1.aapapi[fFsiAch._1.ispin * 3 - 1] * fFsiPrf.aks / 2 - fsi_2cha__.aak;
922 TComplex C10 = GetC(6 + (ISPIN - 1) * 2) + TComplex(a_re, a_im);
923 a_re = fsi_aapap__1.aapapr[fFsiAch._1.ispin * 3 - 1] * AK2S0 / 2;
924 a_im = fsi_aapap__1.aapapi[fFsiAch._1.ispin * 3 - 1] * AK2S0 / 2;
925 TComplex C5 = GetC(6 + (ISPIN - 1) * 2 - 1) + TComplex(a_re, a_im);
926 if (AK2S >= 0) {
927 C5 = C5 - TComplex(0, fsi_2cha__.ak2);
928 } else {
929 C5 = C5 + TComplex(fsi_2cha__.ak2, 0);
930 }
931 C10 = C10 * C5 - GetC(7 + (ISPIN - 1) * 2 - 1) * GetC(7 + (ISPIN - 1) * 2 - 1);
932 SetC(C5 / C10, ISPIN - 1);
933 SetC(-GetC(7 + (ISPIN - 1) * 2 - 1) / C10, ISPIN + 2 - 1);
934 SetC(C5, 5 - 1);
935 SetC(C10, 10 - 1);
936 } /* cpap_ */
937
938 void FemtoFsiParsed::Fsiin(int i_itest__, int i_ich__, int i_iqs__, int i_isi__, int i_i3c__) {
939 /* Initialized data */
940
941 /* Local variables */
942 int j1, j2, itest;
943
944 ++ifirst;
945 /* =======< constants >======================== */
946 fFsiCons.w = 5.0684237202230102;
947 /* from fm to 1/GeV */
948
949 fFsiCons.spi = TMath::Sqrt(TMath::Pi());
950 fFsiCons.dr = 180. / TMath::Pi();
951 /* from radian to degree */
952 fFsiPoc.ac1 = 1e10;
953 fFsiPoc.ac2 = 1e10;
954 /* =======< condition de calculs >============= */
955 // nunit = 11;
956 /* NUNIT=4 ! for SUN in Prague */
957 /* ++ CALL readint4(NUNIT,'ITEST ',ITEST) */
958 /* ++ CALL readint4(NUNIT,'NS ',NS) */
959 /* for IBM in Nantes */
960 itest = i_itest__;
961 if (itest == 1) {
962 /* ++ CALL readint4(NUNIT,'ICH ',ICH) */
963 /* ++ CALL readint4(NUNIT,'IQS ',IQS) */
964 /* ++ CALL readint4(NUNIT,'ISI ',ISI) */
965 /* ++ CALL readint4(NUNIT,'I3C ',I3C) */
966 fFsiNs.ich = i_ich__;
967 fFsiNs.iqs = i_iqs__;
968 fFsiNs.isi = i_isi__;
969 fFsiNs.i3c = i_i3c__;
970 }
971 /* ============================================ */
972 if (ifirst <= 1) {
973 for (j1 = 1; j1 <= 30; ++j1) {
974 for (j2 = 1; j2 <= 10; ++j2) {
975 fsi_fdh__1.fdh[j1 + j2 * 30 - 31] *= fFsiCons.w;
976 fsi_fdh__1.rdh[j1 + j2 * 30 - 31] *= fFsiCons.w;
977 /* L3: */
978 fsi_fdh__1.rbh[j1 + j2 * 30 - 31] *= fFsiCons.w;
979 }
980 }
981 /* print *,"FD,RD,EB,RB: ",FDH(J1,J2),RDH(J1,J2),EBH(J1,J2),RBH(J1,J2)
982 */
983 }
984 /* =================================== */
985 } /* fsiin_ */
986
987 void FemtoFsiParsed::Llini(int lll, int i_ns__, int i_itest__) {
988 /* Initialized data */
989
990 // int ifirst = 0;
991 /* System generated locals */
992 int i__1, i__2, i__3, i__4, i__5, i__6;
993 Double_t d__1, d__2;
994 complex z__1, z__2, z__3, z__4, z__5;
995
996 /* Local variables */
997 Double_t h__;
998 int i3;
999 Double_t c12, hh;
1000 int jh, jj;
1001 Double_t akh;
1002 int msp;
1003 Double_t daks, aksh;
1004 Double_t gpi1h, gpi2h;
1005
1006 /* ===> Initialisation for a given LL value. */
1007 /* =========================================== */
1008 /* ++----- add to be able to call several time------- */
1009 ++ifirst;
1010 /* ===> LL - Initialisation ======================================== */
1011 /* ---- Setting particle masses and charges */
1012 fFsiNs.ll = lll;
1013 fFsiNs.ns = i_ns__;
1014 fFsiPoc.am1 = fsi_amch__1.am1h[fFsiNs.ll - 1];
1015 fFsiPoc.am2 = fsi_amch__1.am2h[fFsiNs.ll - 1];
1016 fFsiPoc.c1 = fsi_amch__1.c1h[fFsiNs.ll - 1];
1017 fFsiPoc.c2 = fsi_amch__1.c2h[fFsiNs.ll - 1];
1018 /* print *,"llini: ",AM1,AM2,C1,C2 */
1019 /* ---> Switches: */
1020 /* ISI=1(0) the strong interaction between the two particles ON (OFF) */
1021 /* IQS=1(0) the quantum statistics ON (OFF); */
1022 /* should be OFF for nonidentical particles */
1023 /* I3C=1(0) the Coulomb interaction with the nucleus ON (OFF) */
1024 /* I3S=1(0) the strong interaction with the nucleus ON (OFF) */
1025 /* ICH=1(0) if C1*C2 is different from 0 (is equal to 0) */
1026 /* - To switch off the Coulomb force between the two particles */
1027 /* put ICH=0 and substitute the strong amplitude parameters by */
1028 /* the ones not affected by Coulomb interaction */
1029 /* ---- --------------------------------------------------------------------
1030 */
1031 if (i_itest__ != 1) {
1032 fFsiNs.ich = 0;
1033 if (fFsiPoc.c1 * fFsiPoc.c2 != 0.) { fFsiNs.ich = 1; }
1034 fFsiNs.iqs = 0;
1035 if (fFsiPoc.c1 + fFsiPoc.am1 == fFsiPoc.c2 + fFsiPoc.am2) { fFsiNs.iqs = 1; }
1036 fFsiNs.i3s = 0;
1037 /* only this option is available */
1038 fFsiNs.isi = 1;
1039 fFsiNs.i3c = 1;
1040 }
1041 /* ---> Calcul. twice the reduced mass (AM), the relative */
1042 /* mass difference (DM) and the Bohr radius (AC) */
1043 fsi_c__._2.am = fFsiPoc.am1 * 2 * fFsiPoc.am2 / (fFsiPoc.am1 + fFsiPoc.am2);
1044 fsi_c__._2.ams = fsi_c__._2.am * fsi_c__._2.am;
1045 fsi_c__._2.dm = (fFsiPoc.am1 - fFsiPoc.am2) / (fFsiPoc.am1 + fFsiPoc.am2);
1046 fFsiAch._1.ac = 1e10;
1047 c12 = fFsiPoc.c1 * fFsiPoc.c2;
1048 if (c12 != 0.) { fFsiAch._1.ac = 274.072 / (c12 * fsi_c__._2.am); }
1049 /* ---Setting spin factors */
1050 fFsiAch._1.mspin = fsi_amch__1.mspinh[fFsiNs.ll - 1];
1051 /* print *,"MSPIN: ",MSPIN */
1052 msp = fFsiAch._1.mspin;
1053 for (fFsiAch._1.ispin = 1; fFsiAch._1.ispin <= 10; ++fFsiAch._1.ispin) {
1054 /* L91: */
1055 fsi_spin__.rho[fFsiAch._1.ispin - 1] = fsi_rhoh__1.rhoh[fFsiNs.ll + fFsiAch._1.ispin * 30 - 31];
1056 }
1057 /* 91 print *,"RHO: ",ISPIN,RHO(ISPIN) */
1058 /* ---> Integration limit AA in the spherical wave approximation */
1059 fsi_aa__.aa = 0.;
1060 /* c IF(NS.EQ.2.OR.NS.EQ.4)AA=.5D0 !!in 1/GeV --> 0.1 fm */
1061 if (fFsiNs.ns == 2 || fFsiNs.ns == 4) { fsi_aa__.aa = 6.; }
1062 /* ---> Setting scatt. length (FD), eff. radius (RD) and, if possible, */
1063 /* -- also the corresp. square well parameters (EB, RB) */
1064 /* !in 1/GeV --> 1.2 fm */
1065 i__1 = msp;
1066 for (jj = 1; jj <= i__1; ++jj) {
1067 fFsiAch._1.ispin = jj;
1068 fsi_fd__.fd[jj - 1] = fsi_fdh__1.fdh[fFsiNs.ll + jj * 30 - 31];
1069 fsi_fd__.rd[jj - 1] = fsi_fdh__1.rdh[fFsiNs.ll + jj * 30 - 31];
1070 fsi_sw__.eb[jj - 1] = fsi_fdh__1.ebh[fFsiNs.ll + jj * 30 - 31];
1071 fsi_sw__.rb[jj - 1] = fsi_fdh__1.rbh[fFsiNs.ll + jj * 30 - 31];
1072 /* print *,"FD,RD,EB,RB: ",FD(JJ),RD(JJ),EB(JJ),RB(JJ) */
1073 /* ---Resets FD and RD for a nucleon-deuteron system (LL=8,9) */
1074 if (fFsiNs.ll == 8 || fFsiNs.ll == 9) {
1075 jh = fFsiNs.ll - 7 + (jj << 1) - 2;
1076 fsi_fd__.fd[jj - 1] = fFsiAand.e_1[jh * 20 - 20];
1077 fsi_fd__.rd[jj - 1] = fFsiAand.e_1[jh * 20 - 19] - fFsiAand.e_1[jh * 20 - 18] * 2 / fFsiAand.e_1[jh * 20 - 20];
1078 }
1079 /* ---Resets FD and RD for a pion-pion system (LL=5,6,7) */
1080 if (fFsiNs.ll == 5 || fFsiNs.ll == 6 || fFsiNs.ll == 7) {
1081 if (fFsiNs.ll == 7) { fsi_fd__.fd[jj - 1] = fFsiAapi.e_1[20] / fsi_c__._2.am; }
1082 if (fFsiNs.ll == 5) { fsi_fd__.fd[jj - 1] = (fFsiAapi.e_1[0] * .6667 + fFsiAapi.e_1[20] * .3333) / fsi_c__._2.am; }
1083 if (fFsiNs.ll == 6) { fsi_fd__.fd[jj - 1] = (fFsiAapi.e_1[0] * .3333 + fFsiAapi.e_1[20] * .6667) / fsi_c__._2.am; }
1084 fFsiPrf.aks = 0.;
1085 daks = 1e-5;
1086 aksh = fFsiPrf.aks + daks;
1087 akh = TMath::Sqrt(aksh);
1088 gpi1h = Gpipi(aksh, c__1);
1089 gpi2h = Gpipi(aksh, c__2);
1090 h__ = 1 / fsi_fd__.fd[jj - 1];
1091 if (fFsiNs.ll == 7) {
1092 i__2 = jj - 1;
1093 d__1 = -akh;
1094 z__2.r = gpi2h, z__2.i = d__1;
1095 z_div(&z__1, &c_b2, &z__2);
1096 fsi_c__._2.c__[i__2].r = z__1.r, fsi_c__._2.c__[i__2].i = z__1.i;
1097 }
1098 if (fFsiNs.ll == 5) {
1099 i__2 = jj - 1;
1100 d__1 = -akh;
1101 z__3.r = gpi1h, z__3.i = d__1;
1102 z_div(&z__2, &c_b33, &z__3);
1103 d__2 = -akh;
1104 z__5.r = gpi2h, z__5.i = d__2;
1105 z_div(&z__4, &c_b34, &z__5);
1106 z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
1107 fsi_c__._2.c__[i__2].r = z__1.r, fsi_c__._2.c__[i__2].i = z__1.i;
1108 }
1109 if (fFsiNs.ll == 6) {
1110 i__2 = jj - 1;
1111 d__1 = -akh;
1112 z__3.r = gpi1h, z__3.i = d__1;
1113 z_div(&z__2, &c_b34, &z__3);
1114 d__2 = -akh;
1115 z__5.r = gpi2h, z__5.i = d__2;
1116 z_div(&z__4, &c_b33, &z__5);
1117 z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
1118 fsi_c__._2.c__[i__2].r = z__1.r, fsi_c__._2.c__[i__2].i = z__1.i;
1119 }
1120 z_div(&z__1, &c_b2, &fsi_c__._2.c__[jj - 1]);
1121 hh = z__1.r;
1122 fsi_fd__.rd[jj - 1] = (hh - h__) * 2 / daks;
1123 }
1124 /* ---Resets FD and RD for a pion-nucleon system (LL=12,13) */
1125 if (fFsiNs.ll == 12 || fFsiNs.ll == 13) {
1126 if (fFsiNs.ll == 12) { fsi_fd__.fd[jj - 1] = fFsiApin.e_1[20]; }
1127 if (fFsiNs.ll == 13) { fsi_fd__.fd[jj - 1] = fFsiApin.e_1[0] * .6667 + fFsiApin.e_1[20] * .3333; }
1128 fFsiPrf.aks = 0.;
1129 daks = 1e-5;
1130 aksh = fFsiPrf.aks + daks;
1131 akh = TMath::Sqrt(aksh);
1132 gpi1h = Gpin(aksh, c__1);
1133 gpi2h = Gpin(aksh, c__2);
1134 h__ = 1 / fsi_fd__.fd[jj - 1];
1135 if (fFsiNs.ll == 12) {
1136 i__2 = jj - 1;
1137 d__1 = -akh;
1138 z__2.r = gpi2h, z__2.i = d__1;
1139 z_div(&z__1, &c_b2, &z__2);
1140 fsi_c__._2.c__[i__2].r = z__1.r, fsi_c__._2.c__[i__2].i = z__1.i;
1141 }
1142 if (fFsiNs.ll == 13) {
1143 i__2 = jj - 1;
1144 d__1 = -akh;
1145 z__3.r = gpi1h, z__3.i = d__1;
1146 z_div(&z__2, &c_b33, &z__3);
1147 d__2 = -akh;
1148 z__5.r = gpi2h, z__5.i = d__2;
1149 z_div(&z__4, &c_b34, &z__5);
1150 z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
1151 fsi_c__._2.c__[i__2].r = z__1.r, fsi_c__._2.c__[i__2].i = z__1.i;
1152 }
1153 z_div(&z__1, &c_b2, &fsi_c__._2.c__[jj - 1]);
1154 hh = z__1.r;
1155 fsi_fd__.rd[jj - 1] = (hh - h__) * 2 / daks;
1156 }
1157 /* ---fm to 1/GeV for pp-bar system */
1158 if (fFsiNs.ll == 30) {
1159 if (ifirst <= 1) {
1160 for (i3 = 1; i3 <= 3; ++i3) {
1161 fsi_aapap__1.aapapr[i3 + jj * 3 - 4] *= fFsiCons.w;
1162 /* L4: */
1163 fsi_aapap__1.aapapi[i3 + jj * 3 - 4] *= fFsiCons.w;
1164 }
1165 /* 4 print *,"AAPAPR,AAPAPI: ",AAPAPR(I3,JJ),AAPAPI(I3,JJ) */
1166 /* --- Calculates complex elements M11=M22=C(6), M12=M21=C(7) for I=0
1167 */
1168 /* --- at k*=0 M11=M22=C(8), M12=M21=C(9) for I=1 */
1169 i__2 = ((jj - 1) << 1) + 6; // suggest parentheses around ‘-’
1170 i__3 = jj * 3 - 3;
1171 i__4 = jj * 3 - 3;
1172 z__3.r = fsi_aapap__1.aapapr[i__3], z__3.i = fsi_aapap__1.aapapi[i__4];
1173 d__1 = 2.;
1174 z__2.r = d__1 * z__3.r, z__2.i = d__1 * z__3.i;
1175 i__5 = jj * 3 - 2;
1176 i__6 = jj * 3 - 2;
1177 z__4.r = fsi_aapap__1.aapapr[i__5], z__4.i = fsi_aapap__1.aapapi[i__6];
1178 z__1.r = z__2.r * z__4.r - z__2.i * z__4.i, z__1.i = z__2.r * z__4.i + z__2.i * z__4.r;
1179 fsi_c__._2.c__[i__2].r = z__1.r, fsi_c__._2.c__[i__2].i = z__1.i;
1180 /* 2a_0Sa_1S */
1181 i__2 = ((jj - 1) << 1) + 5; // suggest parentheses around ‘-’
1182 d__1 = fsi_aapap__1.aapapr[jj * 3 - 3] + fsi_aapap__1.aapapr[jj * 3 - 2];
1183 d__2 = fsi_aapap__1.aapapi[jj * 3 - 3] + fsi_aapap__1.aapapi[jj * 3 - 2];
1184 z__2.r = d__1, z__2.i = d__2;
1185 z_div(&z__1, &z__2,
1186 &fsi_c__._2.c__[((jj - 1) << 1) + 6]); // suggest parentheses around ‘-’
1187 fsi_c__._2.c__[i__2].r = z__1.r, fsi_c__._2.c__[i__2].i = z__1.i;
1188 /* M11=M22 */
1189 i__2 = ((jj - 1) << 1) + 6; // suggest parentheses around ‘-’
1190 d__1 = fsi_aapap__1.aapapr[jj * 3 - 3] - fsi_aapap__1.aapapr[jj * 3 - 2];
1191 d__2 = fsi_aapap__1.aapapi[jj * 3 - 3] - fsi_aapap__1.aapapi[jj * 3 - 2];
1192 z__3.r = d__1, z__3.i = d__2;
1193 z__2.r = -z__3.r, z__2.i = -z__3.i;
1194 z_div(&z__1, &z__2,
1195 &fsi_c__._2.c__[((jj - 1) << 1) + 6]); // suggest parentheses around ‘-’
1196 fsi_c__._2.c__[i__2].r = z__1.r, fsi_c__._2.c__[i__2].i = z__1.i;
1197 /* M12=M21 */
1198 }
1199 }
1200 /* ---Calculation continues for any system (any LL) */
1201 /* L55: */
1202 }
1203 } /* llini_ */
1204
1205 /* ======================================================= */
1206
1207 /* ++ This routine is used to init mass and charge of the nucleus. */
1208 void FemtoFsiParsed::Fsinucl(Double_t r_amn__, Double_t r_cn__) {
1209 fFsiPoc.amn = r_amn__;
1210 fFsiPoc.cn = r_cn__;
1211 } /* fsinucl_ */
1212
1213 /* ====================================================== */
1214
1215 void FemtoFsiParsed::SetMomentum(Double_t* pp1, Double_t* pp2) {
1216 /* particle momenta in NRF */
1217 /* Print *,"momentum",pp1,pp2 */
1218 /* Parameter adjustments */
1219 --pp2;
1220 --pp1;
1221
1222 /* Function Body */
1223 fFsiMom.p1x = pp1[1];
1224 fFsiMom.p1y = pp1[2];
1225 fFsiMom.p1z = pp1[3];
1226 fFsiMom.p2x = pp2[1];
1227 fFsiMom.p2y = pp2[2];
1228 fFsiMom.p2z = pp2[3];
1229 } /* fsimomentum_ */
1230
1231 /* ====================================================== */
1232
1233 void FemtoFsiParsed::SetPosition(Double_t* xt1, Double_t* xt2) {
1234 /* 4-coord. of emis. points in NRF */
1235 /* lc print *,'fsi',xt1,xt2 */
1236 /* Parameter adjustments */
1237 --xt2;
1238 --xt1;
1239
1240 /* Function Body */
1241 fFsiCoor.x1 = xt1[1];
1242 fFsiCoor.y1 = xt1[2];
1243 fFsiCoor.z1 = xt1[3];
1244 fFsiCoor.t1 = xt1[4];
1245 fFsiCoor.x2 = xt2[1];
1246 fFsiCoor.y2 = xt2[2];
1247 fFsiCoor.z2 = xt2[3];
1248 fFsiCoor.t2 = xt2[4];
1249 } /* fsiposition_ */
1250
1251 /* ====================================================== */
1252 /* ====================================================== */
1253
1254 void FemtoFsiParsed::Boosttoprf() {
1255 /* Local variables */
1256 Double_t h1, ts, xs, ys, zs, rs12;
1257
1258 /* part. momenta in NRF */
1259 /* 4-coord. of emis. points in NRF */
1260 xs = fFsiCoor.x1 - fFsiCoor.x2;
1261 ys = fFsiCoor.y1 - fFsiCoor.y2;
1262 zs = fFsiCoor.z1 - fFsiCoor.z2;
1263 ts = fFsiCoor.t1 - fFsiCoor.t2;
1264 rs12 = xs * fSi_p12.p12x + ys * fSi_p12.p12y + zs * fSi_p12.p12z;
1265 h1 = (rs12 / fSi_p12.epm - ts) / fSi_p12.am12;
1266 fFsiPrf.x = xs + fSi_p12.p12x * h1;
1267 fFsiPrf.y = ys + fSi_p12.p12y * h1;
1268 fFsiPrf.z__ = zs + fSi_p12.p12z * h1;
1269 fFsiPrf.t = (fSi_p12.e12 * ts - rs12) / fSi_p12.am12;
1270 fFsiPrf.rps = fFsiPrf.x * fFsiPrf.x + fFsiPrf.y * fFsiPrf.y + fFsiPrf.z__ * fFsiPrf.z__;
1271 fFsiPrf.rp = TMath::Sqrt(fFsiPrf.rps);
1272 /* W WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T */
1273 /* L38: */
1274 fsi_cvk__.cvk =
1275 (fSi_p12.p12x * fFsiPrf.ppx + fSi_p12.p12y * fFsiPrf.ppy + fSi_p12.p12z * fFsiPrf.ppz) / (fSi_p12.p12 * fFsiPrf.ak);
1276 fsi_cvk__.v = fSi_p12.p12 / fSi_p12.e12;
1277 return;
1278 } /* boosttoprf_ */
1279
1280 void FemtoFsiParsed::Fsiwf(Double_t& wei) {
1281 /* System generated locals */
1282 int i__1, i__2, i__3;
1283 Double_t d__1;
1284 complex z__1, z__2, z__3, z__4, z__5;
1285
1286 /* Local variables */
1287 Double_t a1, a2, h1;
1288 int jj, jh;
1289 Double_t xh, ts, xs, ys, zs;
1290 Double_t gak, akf;
1291 Double_t hcp, hra, rs12, xra;
1292 int msp;
1293 Double_t gpi2, gpi1;
1294 Double_t akach;
1295
1296 /* ==> Prepares necessary quantities and call VZ(WEI) to calculate */
1297 /* the weight due to FSI */
1298 /* part. momenta in NRF */
1299 /* 4-coord. of emis. points in NRF */
1300 /* ==>calculating relative 4-coordinates of the particles in PRF */
1301 /* - {T,X,Y,Z} from the relative coordinates {TS,XS,YS,ZS} in NRF */
1302 /* k* (kappa) for 2-nd */
1303 xs = fFsiCoor.x1 - fFsiCoor.x2;
1304 ys = fFsiCoor.y1 - fFsiCoor.y2;
1305 zs = fFsiCoor.z1 - fFsiCoor.z2;
1306 ts = fFsiCoor.t1 - fFsiCoor.t2;
1307 rs12 = xs * fSi_p12.p12x + ys * fSi_p12.p12y + zs * fSi_p12.p12z;
1308 h1 = (rs12 / fSi_p12.epm - ts) / fSi_p12.am12;
1309 fFsiPrf.x = xs + fSi_p12.p12x * h1;
1310 fFsiPrf.y = ys + fSi_p12.p12y * h1;
1311 fFsiPrf.z__ = zs + fSi_p12.p12z * h1;
1312 fFsiPrf.t = (fSi_p12.e12 * ts - rs12) / fSi_p12.am12;
1313 fFsiPrf.rps = fFsiPrf.x * fFsiPrf.x + fFsiPrf.y * fFsiPrf.y + fFsiPrf.z__ * fFsiPrf.z__;
1314 fFsiPrf.rp = TMath::Sqrt(fFsiPrf.rps);
1315 /* W WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T */
1316 /* L38: */
1317 fsi_cvk__.cvk =
1318 (fSi_p12.p12x * fFsiPrf.ppx + fSi_p12.p12y * fFsiPrf.ppy + fSi_p12.p12z * fFsiPrf.ppz) / (fSi_p12.p12 * fFsiPrf.ak);
1319 fsi_cvk__.v = fSi_p12.p12 / fSi_p12.e12;
1320 /* ACH=1.D0 */
1321 if (fFsiNs.ich == 0) { goto L21; }
1322 xh = fFsiAch._1.ac * fFsiPrf.ak;
1323 fFsiAch._1.ach = Acp(xh);
1324 fFsiAch._1.achr = TMath::Sqrt(fFsiAch._1.ach);
1325 fEta = 0.;
1326 if (xh != 0.) { fEta = 1 / xh; }
1327 /* ---HCP, HPR needed (e.g. in GST) if ICH=1 */
1328 hcp = Hc(xh);
1329 fFsiAch._1.hpr = hcp + .1544313298;
1330 /* AAK=ACH*AK ! */
1331 /* HCP2=2*HCP/AC ! needed to calculate C(JJ) for charged particles */
1332 L21:
1333 msp = fFsiAch._1.mspin;
1334 i__1 = msp;
1335 for (jj = 1; jj <= i__1; ++jj) {
1336 fFsiAch._1.ispin = jj;
1337 if (fFsiNs.ns != 1) { goto L22; }
1338 /* ---Calc. quantities for the square well potential; */
1339 /* -- for LL=6-26 the square well potential is not possible or available */
1340 if (fFsiNs.ll == 4) { goto L22; }
1341 /* Computing 2nd power */
1342 d__1 = fsi_sw__.eb[jj - 1];
1343 fsi_sw__.bk[jj - 1] = TMath::Sqrt(d__1 * d__1 + fFsiPrf.aks);
1344 xra = fsi_sw__.rb[jj - 1] * 2 / fFsiAch._1.ac;
1345 hra = fsi_sw__.bk[jj - 1] * fsi_sw__.rb[jj - 1];
1346 Seq(xra, hra);
1347 fsi_sw__.sbkrb[jj - 1] = hra * fsi_bp__.b;
1348 hra = fFsiPrf.ak * fsi_sw__.rb[jj - 1];
1349 Gst(xra, hra);
1350 fsi_sw__.sdk[jj - 1] = fsi_shh__.sh;
1351 fsi_sw__.cdk[jj - 1] = fsi_shh__.chh;
1352 fsi_sw__.sdkk[jj - 1] = fsi_sw__.rb[jj - 1];
1353 if (fFsiPrf.ak != 0.) { fsi_sw__.sdkk[jj - 1] = fsi_shh__.sh / fFsiPrf.ak; }
1354 if (fFsiNs.ich == 1) { fsi_sw__.sdk[jj - 1] = fFsiAch._1.ach * fsi_sw__.sdk[jj - 1]; }
1355 L22:
1356 /* -----------------------------------------------------------------------
1357 */
1358 /* ---Calc. the strong s-wave scattering amplitude = C(JJ) */
1359 /* -- divided by Coulomb penetration factor squared (if ICH=1) */
1360 if (fFsiNs.ns != 1) { goto L230; }
1361 if (fFsiNs.ll != 4) { goto L230; }
1362 /* SW scat. amplitude used for alfa-alfa only */
1363 gak = G(fFsiPrf.ak);
1364 akach = fFsiPrf.ak;
1365 if (fFsiNs.ich == 1) { akach = fFsiPrf.ak * fFsiAch._1.ach; }
1366 i__2 = jj - 1;
1367 d__1 = -akach;
1368 z__2.r = gak, z__2.i = d__1;
1369 z_div(&z__1, &c_b2, &z__2);
1370 fsi_c__._2.c__[i__2].r = z__1.r, fsi_c__._2.c__[i__2].i = z__1.i;
1371 /* amplitude for the SW-potential */
1372 goto L30;
1373 L230:
1374 if (fFsiNs.ll == 5 || fFsiNs.ll == 6 || fFsiNs.ll == 7) { goto L20; }
1375 /* pipi */
1376 if (fFsiNs.ll == 12 || fFsiNs.ll == 13) { goto L20; }
1377 /* piN */
1378 if (fFsiNs.ll == 8 || fFsiNs.ll == 9 || fFsiNs.ll == 18) { goto L20; }
1379 /* Nd, dd */
1380 if (fFsiNs.ll == 14 || fFsiNs.ll == 17 || fFsiNs.ll == 23) { goto L27; }
1381 /* K+K-, K-p, K0K0-b */
1382 if (fFsiNs.ll == 30) { goto L28; }
1383 /* pp-bar */
1384 a1 = fsi_fd__.rd[jj - 1] * fsi_fd__.fd[jj - 1] * fFsiPrf.aks;
1385 a2 = a1 * .5 + 1;
1386 if (fFsiNs.ich == 1) { a2 -= hcp * 2 * fsi_fd__.fd[jj - 1] / fFsiAch._1.ac; }
1387 akf = fFsiPrf.ak * fsi_fd__.fd[jj - 1];
1388 if (fFsiNs.ich == 1) { akf *= fFsiAch._1.ach; }
1389 i__2 = jj - 1;
1390 i__3 = jj - 1;
1391 z__2.r = fsi_fd__.fd[i__3], z__2.i = 0.;
1392 d__1 = -akf;
1393 z__3.r = a2, z__3.i = d__1;
1394 z_div(&z__1, &z__2, &z__3);
1395 fsi_c__._2.c__[i__2].r = z__1.r, fsi_c__._2.c__[i__2].i = z__1.i;
1396 goto L30;
1397 L20:
1398 /* ---Calc. scatt. ampl. C(JJ) for pipi, piN and Nd, dd */
1399 jh = fFsiNs.ll - 7 + (jj << 1) - 2;
1400 if (fFsiNs.ll == 8 || fFsiNs.ll == 9) { gpi2 = Gnd(fFsiPrf.aks, jh); }
1401 if (fFsiNs.ll == 18) { gpi2 = Gdd(fFsiPrf.aks, jj); }
1402 if (fFsiNs.ll == 5 || fFsiNs.ll == 6 || fFsiNs.ll == 7) { gpi2 = Gpipi(fFsiPrf.aks, c__2); }
1403 if (fFsiNs.ll == 12 || fFsiNs.ll == 13) { gpi2 = Gpin(fFsiPrf.aks, c__2); }
1404 i__2 = jj - 1;
1405 d__1 = -fFsiPrf.ak;
1406 z__2.r = gpi2, z__2.i = d__1;
1407 z_div(&z__1, &c_b2, &z__2);
1408
1409 fsi_c__._2.c__[i__2].r = z__1.r,
1410 fsi_c__._2.c__[i__2].i = z__1.i; // z__1 jest ok
1411 /* pi+pi+, nd, pd, pi+p, dd */
1412 if (fFsiNs.ll != 5 && fFsiNs.ll != 6 && fFsiNs.ll != 13) { goto L27; }
1413 if (fFsiNs.ll == 5 || fFsiNs.ll == 6) { gpi1 = Gpipi(fFsiPrf.aks, c__1); }
1414 if (fFsiNs.ll == 13) { gpi1 = Gpin(fFsiPrf.aks, c__1); }
1415 if (fFsiNs.ll == 5 || fFsiNs.ll == 13) {
1416 i__2 = jj - 1;
1417 d__1 = -fFsiPrf.ak;
1418 z__3.r = gpi1, z__3.i = d__1;
1419 z_div(&z__2, &c_b33, &z__3);
1420 i__3 = jj - 1;
1421 z__4.r = fsi_c__._2.c__[i__3].r * .3333, z__4.i = fsi_c__._2.c__[i__3].i * .3333;
1422 z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
1423 fsi_c__._2.c__[i__2].r = z__1.r, fsi_c__._2.c__[i__2].i = z__1.i;
1424 }
1425 /* pi+pi-,pi */
1426 if (fFsiNs.ll == 6) {
1427 i__2 = jj - 1;
1428 d__1 = -fFsiPrf.ak;
1429 z__3.r = gpi1, z__3.i = d__1;
1430 z_div(&z__2, &c_b34, &z__3);
1431 i__3 = jj - 1;
1432 z__4.r = fsi_c__._2.c__[i__3].r * .6667, z__4.i = fsi_c__._2.c__[i__3].i * .6667;
1433 z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
1434 fsi_c__._2.c__[i__2].r = z__1.r, fsi_c__._2.c__[i__2].i = z__1.i;
1435 }
1436 /* pi0pi0 */
1437 L27:
1438 /* ---Calc. K+K-, K0K0-b or K-p s-wave scatt. ampl. */
1439 if (fFsiNs.ll == 14 || fFsiNs.ll == 23) { Ckkb(); }
1440 /* IF(LL.EQ.17)C(JJ)=DCMPLX(-3.29D0,3.55D0) ! Martin'76 (K-p) */
1441 /* IF(LL.EQ.17)C(JJ)=DCMPLX(3.29D0,3.55D0) ! Martin'76 (K-p) WRONG
1442 * SIGN!!! */
1443 if (fFsiNs.ll == 17) {
1444 i__2 = jj - 1;
1445 fsi_c__._2.c__[i__2].r = -2.585, fsi_c__._2.c__[i__2].i = 4.156;
1446 }
1447 /* IF(LL.EQ.17)C(JJ)=DCMPLX(-3.371D0,3.244D0) ! Martin'81 (K-p) */
1448 /* ---Calc. pi+pi-, pi+pi+, pd, pi+p, pi-p, K+K- or K-p s-wave scatt. ampl.
1449 */
1450 /* -- divided by Coulomb penetration factor squared (if ICH=1) */
1451 /* Borasoy'04 (K-p) */
1452 if (fFsiNs.ich == 0) { goto L30; }
1453 fsi_2cha__.aak = fFsiAch._1.ach * fFsiPrf.ak;
1454
1455 fsi_2cha__.hcp2 = hcp * 2 / fFsiAch._1.ac;
1456 /* needed to calculate C(JJ) for charged particles */
1457 i__2 = jj - 1;
1458
1459 z_div(&z__4, &c_b2, &fsi_c__._2.c__[jj - 1]);
1460
1461 z__3.r = z__4.r - fsi_2cha__.hcp2, z__3.i = z__4.i;
1462 d__1 = fFsiPrf.ak - fsi_2cha__.aak;
1463 z__5.r = 0., z__5.i = d__1;
1464 z__2.r = z__3.r + z__5.r, z__2.i = z__3.i + z__5.i;
1465 z_div(&z__1, &c_b2, &z__2);
1466 fsi_c__._2.c__[i__2].r = z__1.r,
1467 fsi_c__._2.c__[i__2].i = z__1.i; // tu jest ok
1468 goto L30;
1469 L28:
1470 /* ---Calc. pp-bar s-wave scatt. ampl. */
1471 Cpap();
1472 L30:;
1473 }
1474 /* *********************************************************************** */
1475 Vz(wei);
1476 return;
1477 } /* fsiwf_ */
1478
1479 void FemtoFsiParsed::Vz(Double_t& wei) // chyba cos nie bundli tutaj
1480 {
1481 /* System generated locals */
1482 int i__1, i__2, i__3;
1483 Double_t r__1;
1484 Double_t d__1, d__2;
1485 complex q__1;
1486 complex z__1, z__2;
1487
1488 /* Local variables */
1489 complex g;
1490 complex z8;
1491 Double_t hh;
1492 int jj;
1493 Double_t hs;
1494 int msp;
1495 complex psi12, psi21;
1496 Double_t rhos;
1497 int jsign;
1498
1499 /* ==> Calculates the weight WEI due to FSI */
1500 /* k* (kappa) for 2-nd */
1501 wei = 0.;
1502 if (fJrat == 1) { goto L11; }
1503 rhos = fFsiPrf.ak * fFsiPrf.rp;
1504 hs = fFsiPrf.x * fFsiPrf.ppx + fFsiPrf.y * fFsiPrf.ppy + fFsiPrf.z__ * fFsiPrf.ppz;
1505 if (rhos < 15. && rhos + TMath::Abs(hs) < 20.) { goto L2; }
1506 /* ---Calc. EIDC=exp(i*Coul.Ph.); */
1507 /* -- used in calc. of hypergeom. f-s in SEQA, FAS at k*R > 15, 20 */
1508 // std::cout<<"fphiA\t"<<fEta<<std::endl;//OK
1509 r__1 = (Double_t) fEta;
1510 q__1.r = (Double_t) 1., q__1.i = r__1;
1511 z8.r = q__1.r, z8.i = q__1.i;
1512
1513 Cgamma(q__1, z8);
1514 z8.r = q__1.r, z8.i = q__1.i;
1515 r__1 = c_abs(&z8);
1516 // std::cout<<"Z8 cabs"<<z8.r<<" "<<z8.i<<" x "<<r__1<<std::endl;
1517 q__1.r = z8.r / r__1, q__1.i = z8.i / r__1;
1518 // std::cout<<"\t"<<fsi_coulph__.eidc.r<<"
1519 //"<<fsi_coulph__.eidc.i<<std::endl;
1520 fsi_coulph__.eidc.r = q__1.r, fsi_coulph__.eidc.i = q__1.i;
1521 L2:
1522 Ff(rhos, hs);
1523 L11:
1524 msp = fFsiAch._1.mspin;
1525 if (fFsiNs.isi == 0) { goto L4; }
1526 /* the strong interaction ON (OFF) if ISI=1(0) */
1527 if (fFsiPrf.rp < fsi_aa__.aa) { goto L4; }
1528 if (fJrat != 1) { Firt(); }
1529 if (fFsiNs.iqs == 0) { goto L5; }
1530 /* the quantum statistics ON (OFF) if IQS=1(0) */
1531 jsign = -1;
1532 i__1 = msp;
1533 for (jj = 1; jj <= i__1; ++jj) {
1534 jsign = -jsign;
1535 i__2 = jj - 1;
1536 i__3 = jj - 1;
1537 z__1.r = fsi_rr__.f[i__2].r * fsi_c__._2.c__[i__3].r - fsi_rr__.f[i__2].i * fsi_c__._2.c__[i__3].i,
1538 z__1.i = fsi_rr__.f[i__2].r * fsi_c__._2.c__[i__3].i + fsi_rr__.f[i__2].i * fsi_c__._2.c__[i__3].r;
1539 g.r = z__1.r, g.i = z__1.i;
1540 if (fFsiNs.ich == 1) {
1541 z__1.r = fFsiAch._1.achr * g.r, z__1.i = fFsiAch._1.achr * g.i;
1542 g.r = z__1.r, g.i = z__1.i;
1543 }
1544 z__2.r = fsi_fff__.f12.r + g.r, z__2.i = fsi_fff__.f12.i + g.i;
1545 z__1.r = fFsiFppn.ff12.r * z__2.r - fFsiFppn.ff12.i * z__2.i, z__1.i = fFsiFppn.ff12.r * z__2.i + fFsiFppn.ff12.i * z__2.r;
1546 psi12.r = z__1.r, psi12.i = z__1.i;
1547 z__2.r = fsi_fff__.f21.r + g.r, z__2.i = fsi_fff__.f21.i + g.i;
1548 z__1.r = fFsiFppn.ff21.r * z__2.r - fFsiFppn.ff21.i * z__2.i, z__1.i = fFsiFppn.ff21.r * z__2.i + fFsiFppn.ff21.i * z__2.r;
1549 psi21.r = z__1.r, psi21.i = z__1.i;
1550 d__1 = (Double_t) jsign;
1551 z__2.r = d__1 * psi21.r, z__2.i = d__1 * psi21.i;
1552 z__1.r = psi12.r + z__2.r, z__1.i = psi12.i + z__2.i;
1553 g.r = z__1.r, g.i = z__1.i;
1554 /* L1: */
1555 /* Computing 2nd power */
1556 d__1 = g.r;
1557 /* Computing 2nd power */
1558 d__2 = d_imag(&g);
1559 wei += fsi_spin__.rho[jj - 1] * (d__1 * d__1 + d__2 * d__2);
1560 }
1561 goto L8;
1562 L5:
1563 i__1 = msp;
1564 for (jj = 1; jj <= i__1; ++jj) {
1565 i__2 = jj - 1;
1566 i__3 = jj - 1;
1567 z__1.r = fsi_rr__.f[i__2].r * fsi_c__._2.c__[i__3].r - fsi_rr__.f[i__2].i * fsi_c__._2.c__[i__3].i,
1568 z__1.i = fsi_rr__.f[i__2].r * fsi_c__._2.c__[i__3].i + fsi_rr__.f[i__2].i * fsi_c__._2.c__[i__3].r;
1569 g.r = z__1.r, g.i = z__1.i;
1570 if (fFsiNs.ich == 1) {
1571 z__1.r = fFsiAch._1.achr * g.r, z__1.i = fFsiAch._1.achr * g.i;
1572 g.r = z__1.r, g.i = z__1.i;
1573 }
1574 /* W WRITE(6,38)'JJ ',JJ,'F ',DREAL(F(JJ)),DIMAG(F(JJ)) */
1575 /* W WRITE(6,38)'JJ ',JJ,'C ',DREAL(C(JJ)),DIMAG(C(JJ)) */
1576 /* W WRITE(6,38)'JJ ',JJ,'G ',DREAL(G),DIMAG(G) */
1577 /* W WRITE(6,38)'JJ ',JJ,'F12+G ',DREAL(F12+G),DIMAG(F12+G) */
1578 /* W WRITE(6,38)'JJ ',JJ,'F21+G ',DREAL(F21+G),DIMAG(F21+G) */
1579 /* L38: */
1580 z__2.r = fsi_fff__.f12.r + g.r, z__2.i = fsi_fff__.f12.i + g.i;
1581 z__1.r = fFsiFppn.ff12.r * z__2.r - fFsiFppn.ff12.i * z__2.i, z__1.i = fFsiFppn.ff12.r * z__2.i + fFsiFppn.ff12.i * z__2.r;
1582 psi12.r = z__1.r, psi12.i = z__1.i;
1583 /* L6: */
1584 /* Computing 2nd power */
1585 d__1 = psi12.r;
1586 /* Computing 2nd power */
1587 d__2 = d_imag(&psi12);
1588 wei += fsi_spin__.rho[jj - 1] * (d__1 * d__1 + d__2 * d__2);
1589 }
1590 /* --- Account for nn-bar->pp-bar channel --------------------------- */
1591 if (fFsiNs.ll == 30) {
1592 i__1 = msp;
1593 for (jj = 1; jj <= i__1; ++jj) {
1594 i__2 = jj + 1;
1595 /* Computing 2nd power */
1596 d__1 = fsi_c__._2.c__[i__2].r;
1597 /* Computing 2nd power */
1598 d__2 = d_imag(&fsi_c__._2.c__[jj + 1]);
1599 hh = fsi_spin__.rho[jj - 1] * (d__1 * d__1 + d__2 * d__2) * fsi_2cha__.amu2_amu1__ * fFsiAch._1.ach / fFsiPrf.rps;
1600 if (fsi_2cha__.ak2s < 0.) { hh *= exp(fFsiPrf.rp * -2 * fsi_2cha__.ak2); }
1601 /* L61: */
1602 wei += hh;
1603 }
1604 }
1605 /* ------------------------------------------------------------------ */
1606 return;
1607 L4:
1608 z__1.r = fFsiFppn.ff12.r * fsi_fff__.f12.r - fFsiFppn.ff12.i * fsi_fff__.f12.i,
1609 z__1.i = fFsiFppn.ff12.r * fsi_fff__.f12.i + fFsiFppn.ff12.i * fsi_fff__.f12.r;
1610 psi12.r = z__1.r, psi12.i = z__1.i;
1611 // std::cout<<"C++"<<fsi_fff__.f12.r<<"+"<<fsi_fff__.f12.i<<std::endl;///źle
1612 if (fFsiNs.iqs == 0) { goto L50; }
1613 /* the quantum statistics ON (OFF) if IQS=1(0) */
1614 z__1.r = fFsiFppn.ff21.r * fsi_fff__.f21.r - fFsiFppn.ff21.i * fsi_fff__.f21.i,
1615 z__1.i = fFsiFppn.ff21.r * fsi_fff__.f21.i + fFsiFppn.ff21.i * fsi_fff__.f21.r;
1616 psi21.r = z__1.r, psi21.i = z__1.i;
1617 jsign = -1;
1618 i__1 = msp;
1619 for (jj = 1; jj <= i__1; ++jj) {
1620 jsign = -jsign;
1621 d__1 = (Double_t) jsign;
1622 z__2.r = d__1 * psi21.r, z__2.i = d__1 * psi21.i;
1623 z__1.r = psi12.r + z__2.r, z__1.i = psi12.i + z__2.i;
1624 g.r = z__1.r, g.i = z__1.i;
1625 /* L3: */
1626 /* Computing 2nd power */
1627 d__1 = g.r;
1628 /* Computing 2nd power */
1629 d__2 = d_imag(&g);
1630 wei += fsi_spin__.rho[jj - 1] * (d__1 * d__1 + d__2 * d__2);
1631 }
1632 goto L8;
1633 L50:
1634 /* Computing 2nd power */
1635 d__1 = psi12.r;
1636 /* Computing 2nd power */
1637 d__2 = d_imag(&psi12);
1638
1639 wei = d__1 * d__1 + d__2 * d__2;
1640 // std::cout << d__1 << " " << d__2 << std::endl;
1641 // std::cout << "||>" << *wei << std::endl;
1642 return;
1643 L8:
1644 wei /= 2;
1645 return;
1646 } /* vz_ */
1647
1648 void FemtoFsiParsed::Firt() {
1649 /* System generated locals */
1650 int i__1, i__2, i__3, i__4;
1651 Double_t d__1;
1652 complex z__1, z__2;
1653
1654 /* Local variables */
1655 Double_t h__;
1656 int i__;
1657 Double_t a1, a2, h1, a12, c12, a21;
1658 int jj;
1659 Double_t hm, s12, hp, fc1, fc2;
1660 complex ch1;
1661 Double_t fs1, fs2, shh, shk, xra;
1662 int msp;
1663 Double_t skr;
1664#define rss ((Double_t*) &fFsiPrf + 9)
1665#define tss ((Double_t*) &fFsiPrf + 8)
1666 Double_t tssa, tssh;
1667 /* ---CALC. THE F(JJ) */
1668 /* -- F(JJ)*C(JJ)= DEVIATION OF THE BETHE-SALPETER AMPL. FROM PLANE WAVE */
1669 msp = fFsiAch._1.mspin;
1670 i__1 = msp;
1671 for (jj = 1; jj <= i__1; ++jj) {
1672 if (jj > 1) { goto L3; }
1673 xra = *rss * 2 / fFsiAch._1.ac;
1674 if (fFsiPrf.ak != 0.) { goto L2; }
1675 shk = 1.;
1676 fsi_shh__.sh = 0.;
1677 shh = fsi_shh__.sh;
1678 fsi_shh__.chh = 1 / *rss;
1679 goto L3;
1680 L2:
1681 h__ = fFsiPrf.ak * *rss;
1682 Gst(xra, h__);
1683 fsi_shh__.sh /= *rss;
1684 fsi_shh__.chh /= *rss;
1685 shh = fsi_shh__.sh;
1686 if (fFsiNs.ich == 1) { shh = fFsiAch._1.ach * fsi_shh__.sh; }
1687 L3:
1688 if (fFsiNs.ns == 2) { goto L1; }
1689 /* ---F= ASYMPTOTIC FORMULA (T= 0 APPROX.); NS= 4 */
1690 L6:
1691 i__2 = jj - 1;
1692 z__1.r = fsi_shh__.chh, z__1.i = shh;
1693 fsi_rr__.f[i__2].r = z__1.r, fsi_rr__.f[i__2].i = z__1.i;
1694 if (fFsiNs.ns != 1) { goto L10; }
1695 /* ---F INSIDE THE SQUARE-WELL (T= 0 APPROX.); NS= 1 */
1696 if (*rss >= fsi_sw__.rb[jj - 1]) { goto L10; }
1697 if (fFsiPrf.ak != 0. && jj == 1) { shk = fsi_shh__.sh / fFsiPrf.ak; }
1698 h__ = fsi_sw__.bk[jj - 1] * *rss;
1699 Gst(xra, h__);
1700 skr = fsi_bp__.b * fsi_sw__.bk[jj - 1];
1701 i__2 = jj - 1;
1702 i__3 = jj - 1;
1703 i__4 = jj - 1;
1704 z__2.r = fsi_sw__.cdk[i__3], z__2.i = fsi_sw__.sdk[i__4];
1705 z__1.r = skr * z__2.r, z__1.i = skr * z__2.i;
1706 fsi_rr__.f[i__2].r = z__1.r, fsi_rr__.f[i__2].i = z__1.i;
1707 d__1 = fsi_sw__.sdkk[jj - 1] * skr - shk * fsi_sw__.sbkrb[jj - 1];
1708 z__2.r = d__1, z__2.i = 0.;
1709
1710 z_div(&z__1, &z__2, &fsi_c__._2.c__[jj - 1]);
1711
1712 ch1.r = z__1.r, ch1.i = z__1.i;
1713 i__2 = jj - 1;
1714 i__3 = jj - 1;
1715 z__2.r = fsi_rr__.f[i__3].r + ch1.r, z__2.i = fsi_rr__.f[i__3].i + ch1.i;
1716 i__4 = jj - 1;
1717 z__1.r = z__2.r / fsi_sw__.sbkrb[i__4], z__1.i = z__2.i / fsi_sw__.sbkrb[i__4];
1718 fsi_rr__.f[i__2].r = z__1.r, fsi_rr__.f[i__2].i = z__1.i;
1719 goto L10;
1720 L1:
1721 /* ---F= ASYMPTOTIC FORMULA (T= 0 NOT REQUIRED); NS= 2 */
1722 if (jj > 1) { goto L8; }
1723 if (*tss == 0.) { goto L6; }
1724 tssa = TMath::Abs(*tss);
1725 if (fsi_c__._2.dm != 0.) { goto L11; }
1726 h__ = fsi_c__._2.am * .5 / tssa;
1727 if (fFsiPrf.ak != 0.) { goto L4; }
1728 hm = h__ * fFsiPrf.rps;
1729 if (hm >= 3e15) { goto L6; }
1730 fs1 = DfrSin(hm);
1731 fc1 = DfrCos(hm);
1732 fc2 = fc1;
1733 fs2 = fs1;
1734 goto L5;
1735 L4:
1736 h1 = fFsiPrf.ak * tssa / fsi_c__._2.am;
1737 /* Computing 2nd power */
1738 d__1 = *rss - h1;
1739 hm = h__ * (d__1 * d__1);
1740 /* Computing 2nd power */
1741 d__1 = *rss + h1;
1742 hp = h__ * (d__1 * d__1);
1743 if (hp >= 3e15) { goto L6; }
1744 fs1 = DfrSin(hm);
1745 fc1 = DfrCos(hm);
1746 fs2 = DfrSin(hp);
1747 fc2 = DfrCos(hp);
1748 goto L5;
1749 L11:
1750 fs1 = 0.;
1751 fs2 = 0.;
1752 fc1 = 0.;
1753 fc2 = 0.;
1754 for (i__ = 1; i__ <= 2; ++i__) {
1755 if (i__ == 1) { tssh = tssa * (fsi_c__._2.dm + 1); }
1756 if (i__ == 2) { tssh = tssa * (1 - fsi_c__._2.dm); }
1757 h__ = fsi_c__._2.am * .5 / tssh;
1758 if (fFsiPrf.ak != 0.) { goto L12; }
1759 hm = h__ * fFsiPrf.rps;
1760 if (hm >= 3e15) { goto L6; }
1761 fs1 += DfrSin(hm) / 2;
1762 fc1 += DfrCos(hm) / 2;
1763 if (i__ == 1) { goto L13; }
1764 fc2 = fc1;
1765 fs2 = fs1;
1766 goto L13;
1767 L12:
1768 h1 = fFsiPrf.ak * tssh / fsi_c__._2.am;
1769 /* Computing 2nd power */
1770 d__1 = *rss - h1;
1771 hm = h__ * (d__1 * d__1);
1772 /* Computing 2nd power */
1773 d__1 = *rss + h1;
1774 hp = h__ * (d__1 * d__1);
1775 if (hp >= 3e15) { goto L6; }
1776 fs1 += DfrSin(hm) / 2;
1777 fc1 += DfrCos(hm) / 2;
1778 fs2 += DfrSin(hp) / 2;
1779 fc2 += DfrCos(hp) / 2;
1780 L13:;
1781 }
1782 L5:
1783 c12 = fc1 + fs2;
1784 s12 = fs1 + fc2;
1785 a12 = fs1 - fc2;
1786 a21 = fs2 - fc1;
1787 a2 = (fsi_shh__.chh * (a12 + a21) + fsi_shh__.sh * (a12 - a21)) * .5 + shh;
1788 a1 = (fsi_shh__.chh * (c12 + s12) + fsi_shh__.sh * (c12 - s12)) * .5;
1789 i__2 = jj - 1;
1790 z__2.r = a1, z__2.i = a2;
1791 z__1.r = z__2.r * .3989422, z__1.i = z__2.i * .3989422;
1792 fsi_rr__.f[i__2].r = z__1.r, fsi_rr__.f[i__2].i = z__1.i;
1793 goto L10;
1794 L8:
1795 i__2 = jj - 1;
1796 fsi_rr__.f[i__2].r = fsi_rr__.f[0].r, fsi_rr__.f[i__2].i = fsi_rr__.f[0].i;
1797 L10:;
1798 }
1799 return;
1800 } /* firt_ */
1801 Double_t FemtoFsiParsed::Exf(Double_t x) {
1802 if (x < -15.) { return 0; }
1803 return TMath::Exp(x);
1804 } /* exf_ */
1805
1806 void FemtoFsiParsed::Seq(Double_t x, Double_t h__) {
1807 /* Initialized data */
1808 Double_t BH[3] = {1, x * 0.5, 0};
1809 Double_t PH[3] = {1, 0, 0};
1810 Double_t B = fsi_bp__.b; // TODO
1811 B = 1 + BH[1];
1812 Double_t P = fsi_bp__.p; // TODO
1813 P = 1;
1814 Double_t HS = h__ * h__;
1815 Int_t J = 0;
1816 Bool_t calc = kTRUE;
1817 Double_t Z;
1818
1819 do {
1820 J = J + 1;
1821 BH[2] = (x * BH[1] - HS * BH[0]) / ((J + 1) * (J + 2));
1822 PH[2] = (x * PH[1] - HS * PH[0] - (2 * J + 1) * x * BH[1]) / (J * (J + 1));
1823 B = B + BH[2];
1824 P = P + PH[2];
1825 Z = TMath::Abs(BH[1]) + TMath::Abs(BH[2]) + TMath::Abs(PH[1]) + TMath::Abs(PH[2]);
1826 // std::cout<<"Z"<<Z<<std::endl;
1827 if (Z < 1e-7) {
1828 calc = kFALSE;
1829 fsi_bp__.b = B;
1830 fsi_bp__.p = P;
1831 return;
1832 }
1833 BH[0] = BH[1];
1834 BH[1] = BH[2];
1835 PH[0] = PH[1];
1836 PH[1] = PH[2];
1837 } while (calc);
1838 /*
1839 fsi_bp__.b += BH[2];
1840 fsi_bp__.p += PH[2];
1841
1842 Double_t err = 1e-7;
1843
1844 int j;
1845 Double_t z__, bh[3], ph[3], hs;
1846
1847 // ---CALC. FUNCTIONS B, P (EQS. (17) OF G-K-L-L);
1848 // -- NEEDED TO CALC. THE CONFLUENT HYPERGEOMETRIC FUNCTION GST.
1849 bh[0] = 1.;
1850 ph[0] = 1.;
1851 ph[1] = 0.;
1852 bh[1] = x * .5;
1853 fsi_bp__.b = bh[1] + 1;
1854 fsi_bp__.p = 1.;
1855 hs = h__ * h__;
1856 j = 0;
1857 L2:
1858 ++j;
1859 bh[2] = (x * bh[1] - hs * bh[0]) / ((j + 1) * (j + 2));
1860 ph[2] = (x * ph[1] - hs * ph[0] - ((j << 1) + 1) * x * bh[1]) / (j * (j
1861 + 1));
1862 fsi_bp__.b += bh[2];
1863 fsi_bp__.p += ph[2];
1864 z__ = TMath::Abs(bh[1]) + TMath::Abs(bh[2]) + TMath::Abs(ph[1]) +
1865 TMath::Abs(ph[2]); if (z__ < err) { std::cout<<"SEQ"<<fsi_bp__.b-B<<"
1866 "<<fsi_bp__.p-P<<std::endl; return ;
1867 }
1868 bh[0] = bh[1];
1869 bh[1] = bh[2];
1870 ph[0] = ph[1];
1871 ph[1] = ph[2];
1872 goto L2;
1873 */
1874 } /* seq_ */
1875
1876 void FemtoFsiParsed::Seqa(Double_t h__) {
1877 Double_t ARG = h__ - fEta * TMath::Log(2.0 * h__);
1878 TComplex GST(TMath::Cos(ARG), TMath::Sin(ARG));
1879 TComplex EIDC(fsi_coulph__.eidc.r, fsi_coulph__.eidc.i); // TODO
1880 Double_t ACH = fFsiAch._1.ach;
1881 Double_t ACHR = fFsiAch._1.achr;
1882 GST = ACHR * EIDC * GST;
1883 Double_t CHH = GST.Re();
1884 Double_t SH = GST.Im() / ACH;
1885 Double_t B = SH / h__;
1886 fsi_shh__.chh = CHH;
1887 fsi_shh__.sh = SH;
1888 fsi_bp__.b = B;
1889 //=====================================
1890 /* System generated locals */
1891 Double_t d__1, d__2;
1892 complex z__1, z__2;
1893
1894 /* Local variables */
1895 Double_t arg;
1896 complex gst;
1897
1898 /* ---CALC. FUNCTIONS CHH=REAL(GST), SH=IMAG(GST)/ACH, B=SH/H */
1899 /* -- IN THE ASYMPTOTIC REGION H=K*R >> 1. */
1900 arg = h__ - fEta * TMath::Log(h__ * 2);
1901 d__1 = TMath::Cos(arg);
1902 d__2 = TMath::Sin(arg);
1903 z__1.r = d__1, z__1.i = d__2;
1904 gst.r = z__1.r, gst.i = z__1.i;
1905 z__2.r = fFsiAch._1.achr * fsi_coulph__.eidc.r, z__2.i = fFsiAch._1.achr * fsi_coulph__.eidc.i;
1906 z__1.r = z__2.r * gst.r - z__2.i * gst.i, z__1.i = z__2.r * gst.i + z__2.i * gst.r;
1907 gst.r = z__1.r, gst.i = z__1.i;
1908 fsi_shh__.chh = gst.r;
1909 fsi_shh__.sh = d_imag(&gst) / fFsiAch._1.ach;
1910 fsi_bp__.b = fsi_shh__.sh / h__;
1911 // std::cout<<"SEQA"<< fsi_shh__.chh-CHH<<" "<< fsi_shh__.sh-SH<<"
1912 // "<<fsi_bp__.b - B<<std::endl;
1913 return;
1914 } /* seqa_ */
1915
1916 void FemtoFsiParsed::Ff(Double_t rho, Double_t h__) {
1917 /* System generated locals */
1918 Double_t d__1;
1919 complex z__1, z__2;
1920
1921 /* Local variables */
1922 Double_t c__, s;
1923 // Double_t rhom, rhop;
1924
1925 /* ---Calc. F12, F21; */
1926 /* -- F12= FF0* plane wave, FF0=F*ACHR, */
1927 /* ---F is the confluent hypergeometric function, */
1928 /* -- ACHR=sqrt(ACH), where ACH is the Coulomb factor */
1929 c__ = TMath::Cos(h__);
1930 s = TMath::Sin(h__);
1931 d__1 = -s;
1932 z__1.r = c__, z__1.i = d__1;
1933 fsi_fff__.f12.r = z__1.r, fsi_fff__.f12.i = z__1.i;
1934 z__1.r = c__, z__1.i = s;
1935 fsi_fff__.f21.r = z__1.r, fsi_fff__.f21.i = z__1.i;
1936 // std::cout<<"C++FFF"<<fsi_fff__.f12.r <<" "<<fsi_fff__.f12.i<<std::endl;
1937 // std::cout<<"C++FFF"<<fsi_fff__.f21.r <<" "<<fsi_fff__.f21.i<<std::endl;
1938 if (fFsiNs.ich == 0) { return; }
1939 // rhop = rho + h__;
1940 // rhom = rho - h__;
1941 z__2.r = fsi_fff__.f12.r, z__2.i = fsi_fff__.f12.i;
1942 // std::cout<<"C++FFA >"<<z__2.r<<" "<<z__2.i<<") "<<*rho<<"
1943 // "<<*h__<<std::endl;
1944 Ff0(z__2, rho, h__);
1945 // std::cout<<"C++ out1 < ("<<z__2.r<<" "<<z__2.i<<")"<<std::endl;
1946 z__1.r = z__2.r * fsi_fff__.f12.r - z__2.i * fsi_fff__.f12.i, z__1.i = z__2.r * fsi_fff__.f12.i + z__2.i * fsi_fff__.f12.r;
1947 fsi_fff__.f12.r = z__1.r, fsi_fff__.f12.i = z__1.i;
1948 d__1 = -(h__);
1949 z__2.r = fsi_fff__.f21.r, z__2.i = fsi_fff__.f21.i;
1950 // std::cout<<"C++FF0 <("<<z__1.r<<" "<<z__1.i<<") "<<std::endl;
1951 // std::cout<<"C++FFB >"<<z__2.r<<" "<< z__2.i<<") "<<*rho<<"
1952 // "<<-*h__<<std::endl;
1953 Ff0(z__2, rho, d__1);
1954 // std::cout<<"C++ out2 < ("<<z__2.r<<" "<<z__2.i<<")"<<std::endl;
1955 z__1.r = z__2.r * fsi_fff__.f21.r - z__2.i * fsi_fff__.f21.i, z__1.i = z__2.r * fsi_fff__.f21.i + z__2.i * fsi_fff__.f21.r;
1956 // std::cout<<"C++FF0 <("<<z__1.r<<" "<<z__1.i<<") "<<std::endl;
1957 fsi_fff__.f21.r = z__1.r, fsi_fff__.f21.i = z__1.i;
1958 fsi_fff__.f21.r = z__1.r, fsi_fff__.f21.i = z__1.i;
1959 // std::cout<<"C++ FFA < ("<<fsi_fff__.f12.r<<"
1960 // "<<fsi_fff__.f12.i<<")"<<std::endl; std::cout<<"C++ FFB <
1961 // ("<<fsi_fff__.f21.r<<" "<<fsi_fff__.f21.i<<")"<<std::endl;
1962 return;
1963 } /* ff_ */
1964
1965 void FemtoFsiParsed::Fas(complex& ret_val, Double_t rks) {
1966 TComplex FAS;
1967 Double_t D1, D2;
1968 D1 = TMath::Log(rks) * fEta;
1969 D2 = fEta * fEta / rks;
1970 TComplex ZZ1(TMath::Cos(D1), TMath::Sin(D1));
1971 TComplex EIDC(fsi_coulph__.eidc.r, fsi_coulph__.eidc.i); // TODO
1972 ZZ1 = ZZ1 / EIDC;
1973 FAS = (TComplex(1, -D2)) * ZZ1;
1974 FAS = FAS - TComplex(TMath::Cos(rks), TMath::Sin(rks)) * fEta / rks / ZZ1;
1975 ret_val.r = FAS.Re();
1976 ret_val.i = FAS.Im();
1977 } /* fas_ */
1978
1979 /* Double Complex */ void FemtoFsiParsed::Ff0(complex& ret_val, Double_t rho, Double_t h__) {
1980 /* Initialized data */
1981 // std::cout<<"C+ ff0 in ("<<ret_val->r<<"
1982 // "<<ret_val->i<<")\t"<<*rho<<std::endl;
1983 Double_t err = 1e-5;
1984
1985 /* System generated locals */
1986 int i__1;
1987 Double_t d__1;
1988 complex z__1, z__2;
1989
1990 /* Local variables */
1991 complex a;
1992 int j = 0;
1993 complex s, z__;
1994 Double_t zi, zr;
1995 complex alf;
1996 complex alf1;
1997 Double_t rhop;
1998
1999 /* -- FF0=F*ACHR */
2000 /* -- F is the confluent hypergeometric function */
2001 /* -- (Eq. (15) of G-K-L-L), F= 1 at r* << AC */
2002 /* -- ACHR=sqrt(ACH), where ACH is the Coulomb factor */
2003 s.r = 1., s.i = 0.;
2004 ret_val.r = s.r, ret_val.i = s.i;
2005 rhop = rho + h__;
2006 /* C GOTO 5 ! rejects the approx. calcul. of hyperg. f-ion F */
2007 if (rhop < 20.) { goto L5; }
2008 Fas(z__1, rhop);
2009 ret_val.r = z__1.r, ret_val.i = z__1.i;
2010 /* approx. calc. */
2011 return;
2012 L5:
2013 d__1 = -fEta;
2014 z__1.r = 0., z__1.i = d__1;
2015 alf.r = z__1.r, alf.i = z__1.i;
2016 z__1.r = alf.r - 1, z__1.i = alf.i;
2017 alf1.r = z__1.r, alf1.i = z__1.i;
2018 z__1.r = 0., z__1.i = rhop;
2019 z__.r = z__1.r, z__.i = z__1.i;
2020 j = 0;
2021 L3:
2022 ++j;
2023 z__2.r = alf1.r + j, z__2.i = alf1.i;
2024 i__1 = j * j;
2025 d__1 = (Double_t) i__1;
2026 z__1.r = z__2.r / d__1, z__1.i = z__2.i / d__1;
2027 a.r = z__1.r, a.i = z__1.i;
2028 z__2.r = s.r * a.r - s.i * a.i, z__2.i = s.r * a.i + s.i * a.r;
2029 z__1.r = z__2.r * z__.r - z__2.i * z__.i, z__1.i = z__2.r * z__.i + z__2.i * z__.r;
2030 s.r = z__1.r, s.i = z__1.i;
2031 z__1.r = ret_val.r + s.r, z__1.i = ret_val.i + s.i;
2032 ret_val.r = z__1.r, ret_val.i = z__1.i;
2033 zr = (d__1 = s.r, TMath::Abs(d__1));
2034 zi = (d__1 = d_imag(&s), TMath::Abs(d__1));
2035 if (zr + zi > err) { goto L3; }
2036 z__1.r = fFsiAch._1.achr * ret_val.r, z__1.i = fFsiAch._1.achr * ret_val.i;
2037 ret_val.r = z__1.r, ret_val.i = z__1.i;
2038 // std::cout<<"C+ ff0 out ("<<ret_val->r<<" "<<ret_val->i<<")\t"<<std::endl;
2039 return;
2040 } /* ff0_ */
2041
2042 Double_t FemtoFsiParsed::Hc(Double_t xa) {
2043 /* Initialized data */
2044
2045 Double_t bn[15] = {.08333333333,
2046 .008333333333,
2047 .00396825396825,
2048 .004166666667,
2049 .007575757576,
2050 .02109279609,
2051 .08333333333,
2052 .4432598039,
2053 3.05395433,
2054 26.45621212,
2055 281.4601449,
2056 3607.510546,
2057 54827.58333,
2058 974936.8235,
2059 20052695.8};
2060
2061 /* System generated locals */
2062 int i__1;
2063 Double_t ret_val, d__1;
2064
2065 /* Local variables */
2066 int i__, n;
2067 Double_t s, x, x2, ds, xp;
2068 int ima;
2069
2070 /* ---HC = h-function of Landau-Lifshitz: h(x)=Re[psi(1-i/x)]+ln(x) */
2071 /* -- psi(x) is the digamma function (the logarithmic derivative of */
2072 /* -- the gamma function) */
2073 x = TMath::Abs(xa);
2074 if (x < .33) {
2075 x2 = x * x;
2076 xp = x2;
2077 ret_val = 0.;
2078 ima = 9;
2079 if (x < .1) { ima = 3; }
2080 i__1 = ima;
2081 for (i__ = 1; i__ <= i__1; ++i__) {
2082 ret_val += xp * bn[i__ - 1];
2083 /* L4: */
2084 xp = x2 * xp;
2085 }
2086 return ret_val;
2087 }
2088 /* C IF(X.GE.3.5D0) GO TO 2 */
2089 s = 0.;
2090 n = 0;
2091
2092 do {
2093 ++n;
2094 /* Computing 2nd power */
2095 d__1 = n * x;
2096 ds = 1. / n / (d__1 * d__1 + 1);
2097 s += ds;
2098 } while (ds > 1e-13);
2099 ret_val = s - .5772156649 + TMath::Log(x);
2100 return ret_val;
2101
2102 } /* hc_ */
2103
2104 Double_t FemtoFsiParsed::Acp(Double_t x) {
2105 /* --- ACP = COULOMB PENETRATION FACTOR */
2106 if (x < .05 && x >= 0.) { return 1e-6; }
2107 Double_t y = 6.2831853 / x;
2108 return y / (Exf(y) - 1);
2109 } /* acp_ */
2110
2111 void FemtoFsiParsed::Gst(Double_t x, Double_t h__) {
2112 if (fFsiNs.ich == 1) { // goto 2
2113 if (h__ > 15.0) { // goto 4
2114 Seqa(h__);
2115 } else {
2116 Seq(x, h__);
2117 /* exact calculation */
2118 fsi_shh__.sh = h__ * fsi_bp__.b;
2119 fsi_shh__.chh = fsi_bp__.p + fsi_bp__.b * x * (TMath::Log((TMath::Abs(x))) + fFsiAch._1.hpr);
2120 }
2121 } else {
2122 fsi_shh__.sh = TMath::Sin(h__);
2123 fsi_shh__.chh = TMath::Cos(h__);
2124 fsi_bp__.b = 1.;
2125 if (h__ != 0.) { fsi_bp__.b = fsi_shh__.sh / h__; }
2126 fsi_bp__.p = fsi_shh__.chh;
2127 }
2128
2129 /*
2130
2131 if (fFsiNs.ich == 1) {
2132 goto L2;
2133 }
2134
2135 fsi_shh__.sh = TMath::Sin(h__);
2136 fsi_shh__.chh = TMath::Cos(h__);
2137 fsi_bp__.b = 1.;
2138 if (h__ != 0.) {
2139 fsi_bp__.b = fsi_shh__.sh / h__;
2140 }
2141 fsi_bp__.p = fsi_shh__.chh;
2142 return;
2143 L2:
2144 if (h__ > 15.) {
2145 goto L4;
2146 }
2147
2148 Seq(x, h__);
2149
2150 fsi_shh__.sh = h__ * fsi_bp__.b;
2151 fsi_shh__.chh = fsi_bp__.p + fsi_bp__.b * x * (TMath::Log((TMath::Abs(x))) +
2152 fFsiAch._1.hpr);
2153 return;
2154 L4:
2155 Seqa(x, h__);
2156 return;
2157 */
2158 } /* gst_ */
2159
2160 void FemtoFsiParsed::Ff1(complex& ret_val, Double_t rho, Double_t h__) {
2161 /* System generated locals */
2162 Double_t r__1;
2163 complex q__1;
2164 complex z__1;
2165
2166 /* Local variables */
2167 complex z8;
2168
2169 /* ---FF1=FF0; used for particle-nucleus system */
2170 /* -- FF0=F12*ACHR */
2171 /* -- F12 is the confluent hypergeometric function */
2172 /* -- (Eq. (15) of G-K-L-L), F12= 1 at r* << AC */
2173 /* -- ACHR=sqrt(ACH), where ACH is the Coulomb factor */
2174 ret_val.r = 1., ret_val.i = 0.;
2175 if (fIch == 0) { goto L2; }
2176 if (rho < 15. && rho + h__ < 20.) { goto L2; }
2177 /* ---Calc. EIDC=exp(i*Coul.Ph.); */
2178 /* -- used in calc. of hypergeom. f-s in SEQA, FAS at k*R > 15, 20 */
2179 r__1 = (Double_t) fEta;
2180 q__1.r = (Double_t) 1., q__1.i = r__1;
2181 z8.r = q__1.r, z8.i = q__1.i;
2182 Cgamma(q__1, z8);
2183 // std::cout<<"Z8 cabs"<<z8.r<<" "<<z8.i<<" x "<<q__1.r<<"
2184 // "<<q__1.i<<std::endl;
2185 z8.r = q__1.r, z8.i = q__1.i;
2186 r__1 = c_abs(&z8);
2187 q__1.r = z8.r / r__1, q__1.i = z8.i / r__1;
2188 fsi_coulph__.eidc.r = q__1.r, fsi_coulph__.eidc.i = q__1.i;
2189 L2:
2190 Ff0(z__1, rho, h__);
2191 ret_val.r = z__1.r, ret_val.i = z__1.i;
2192 return;
2193 } /* ff1_ */
2194
2195 Double_t FemtoFsiParsed::G(Double_t ak) {
2196 /* System generated locals */
2197 Double_t ret_val, d__1;
2198
2199 /* Local variables */
2200 Double_t h__, x;
2201 Double_t gg, xh, zz;
2202 Double_t hcp, aks;
2203 Double_t brho;
2204 Double_t brhop, brhos, prhos;
2205
2206 /* ---Used to calculate SW scattering amplitude for alpa-alpha system */
2207 /* -- and for sear.f (square well potential search) */
2208 /* ---NOTE THAT SCATT. AMPL.= 1/CMPLX(G(AK),-AK*ACH) */
2209 fsi_dak__.hac = 0.;
2210 fFsiAch._2.ach = 1.;
2211 if (fFsiNs.ich != 0) {
2212 xh = fFsiAch._2.ac * ak;
2213 hcp = Hc(xh);
2214 fFsiAch._2.hpr = hcp + .1544313298;
2215 fFsiAch._2.ach = Acp(xh);
2216 fsi_dak__.hac = hcp * 2 / fFsiAch._2.ac;
2217 }
2218 /* Computing 2nd power */
2219 d__1 = ak;
2220 aks = d__1 * d__1;
2221 /* Computing 2nd power */
2222 d__1 = fsi_sw__.eb[fFsiAch._2.jj - 1];
2223 fsi_sw__.bk[fFsiAch._2.jj - 1] = TMath::Sqrt(aks + d__1 * d__1);
2224 /* kappa=kp */
2225 x = fsi_sw__.rb[fFsiAch._2.jj - 1] * 2 / fFsiAch._2.ac;
2226 h__ = fsi_sw__.bk[fFsiAch._2.jj - 1] * fsi_sw__.rb[fFsiAch._2.jj - 1];
2227 /* kp*d */
2228 Gst(x, h__);
2229 brho = fsi_bp__.b;
2230 /* B(kp,d) */
2231 fsi_sw__.sbkrb[fFsiAch._2.jj - 1] = fsi_shh__.sh;
2232 /* kp*d*B(kp,d) */
2233 Deriw(x, h__);
2234 brhop = fsi_deriv__.bpr;
2235 /* B'(kp,d)= dB(kp,r)/dln(r) at r=d */
2236 h__ = ak * fsi_sw__.rb[fFsiAch._2.jj - 1];
2237 Gst(x, h__);
2238 fsi_sw__.cdk[fFsiAch._2.jj - 1] = fsi_shh__.chh;
2239 /* ReG(k,d) */
2240 brhos = fsi_bp__.b;
2241 /* B(k,d) */
2242 prhos = fsi_bp__.p;
2243 /* P(k,d) */
2244 fsi_sw__.sdk[fFsiAch._2.jj - 1] = fsi_shh__.sh;
2245 if (fFsiNs.ich != 0) {
2246 fsi_sw__.sdk[fFsiAch._2.jj - 1] = fFsiAch._2.ach * fsi_shh__.sh;
2247 /* ImG(k,d) */
2248 if (ak == 0. && fFsiAch._2.ac < 0.) { fsi_sw__.sdk[fFsiAch._2.jj - 1] = x * (Double_t) 3.14159 * fsi_bp__.b; }
2249 }
2250 fsi_sw__.sdkk[fFsiAch._2.jj - 1] = fsi_sw__.rb[fFsiAch._2.jj - 1];
2251 if (ak != 0.) { fsi_sw__.sdkk[fFsiAch._2.jj - 1] = fsi_shh__.sh / ak; }
2252 /* d*B(k,d) */
2253 Deriw(x, h__);
2254 /* PPR=P'(k,d)= dP(k,r)/dln(r) at r=d */
2255 zz = fsi_deriv__.ppr - prhos;
2256 if (fFsiNs.ich == 1) { zz += x * (brhos + fsi_deriv__.bpr * (TMath::Log((TMath::Abs(x))) + fFsiAch._2.hpr)); }
2257 /* ZZ= P'(k,d)-P(k,d)+x*{B(k,d)+B'(k,d)*[ln!x!+2*C-1+h(k*ac)]} */
2258 gg = (brhop * fsi_sw__.cdk[fFsiAch._2.jj - 1] - brho * zz) / fsi_sw__.rb[fFsiAch._2.jj - 1];
2259 /* GG= [B'(kp,d)*ReG(k,d)-B(kp,d)*ZZ]/d */
2260 ret_val = gg / (brho * fsi_deriv__.bpr - brhop * brhos);
2261 /* G= GG/[B(kp,d)*B'(k,d)-B'(kp,d)*B(k,d)] */
2262 return ret_val;
2263 } /* g_ */
2264
2265 void FemtoFsiParsed::Deriw(Double_t x, Double_t h__) {
2266 /* System generated locals */
2267 Double_t d__1;
2268
2269 /* Local variables */
2270 Double_t b1, q1, hh, hhh;
2271
2272 /* ---CALLED BY F-N G(AK) */
2273 hh = 1e-4;
2274 d__1 = h__ - hh;
2275 Gst(x, d__1);
2276 q1 = fsi_bp__.p;
2277 b1 = fsi_bp__.b;
2278 d__1 = h__ + hh;
2279 Gst(x, d__1);
2280 hhh = hh + hh;
2281 fsi_deriv__.bpr = h__ * (fsi_bp__.b - b1) / hhh;
2282 fsi_deriv__.ppr = h__ * (fsi_bp__.p - q1) / hhh;
2283 if (fFsiNs.ich == 0) return;
2284 d__1 = x - hh;
2285 Gst(d__1, h__);
2286 q1 = fsi_bp__.p;
2287 b1 = fsi_bp__.b;
2288 d__1 = x + hh;
2289 Gst(d__1, h__);
2290 fsi_deriv__.bpr += x * (fsi_bp__.b - b1) / hhh;
2291 fsi_deriv__.ppr += x * (fsi_bp__.p - q1) / hhh;
2292 } /* deriw_ */
2293
2294 /* ================================================================ */
2295 /*
2296 int FemtoFsiParsed::read_file__(char *key, char *ch8, int *int4,
2297 Double_t *real8, int *ierr, int *nunit, int key_len,
2298 int ch8_len)
2299 {
2300 //std::cout<<"READING"<<std::endl;
2301 //
2302 // System generated locals
2303 alist al__1, al__2;
2304
2305
2306 Local variables
2307 static char test[10], type__[4];
2308
2309 // Fortran I/O blocks
2310 static cilist io___167 = { 0, 0, 0, "(A10,2X,A4)", 0 };
2311 static cilist io___170 = { 0, 0, 0, "(18X,A8,54X)", 0 };
2312 static cilist io___171 = { 0, 0, 0, "(18X,I8,54X)", 0 };
2313 static cilist io___172 = { 0, 0, 0, "(18X,F10.5,52X)", 0 };
2314
2315
2316 *ierr = 0;
2317 al__1.aerr = 0;
2318 al__1.aunit = *nunit;
2319 f_rew(&al__1);
2320 L1:
2321 io___167.ciunit = *nunit;
2322 s_rsfe(&io___167);
2323 do_fio(&c__1, test, (int)10);
2324 do_fio(&c__1, type__, (int)4);
2325 e_rsfe();
2326 if (s_cmp(test, key, (int)10, (int)10) == 0) {
2327 al__2.aerr = 0;
2328 al__2.aunit = *nunit;
2329 f_back(&al__2);
2330 if (s_cmp(type__, "CHAR", (int)4, (int)4) == 0) {
2331 io___170.ciunit = *nunit;
2332 s_rsfe(&io___170);
2333 do_fio(&c__1, ch8, (int)8);
2334 e_rsfe();
2335 }
2336 if (s_cmp(type__, "INT4", (int)4, (int)4) == 0) {
2337 io___171.ciunit = *nunit;
2338 s_rsfe(&io___171);
2339 do_fio(&c__1, (char *)&(*int4), (int)sizeof(int));
2340 e_rsfe();
2341 }
2342 if (s_cmp(type__, "REA8", (int)4, (int)4) == 0) {
2343 io___172.ciunit = *nunit;
2344 s_rsfe(&io___172);
2345 do_fio(&c__1, (char *)&(*real8), (int)sizeof(Double_t));
2346 e_rsfe();
2347 }
2348 } else {
2349 if (s_cmp(test, "* E.O.F. *", (int)10, (int)10) != 0) {
2350 goto L1;
2351 } else {
2352 *ierr = 1;
2353 }
2354 }
2355
2356 // IF(IERR.EQ.1)STOP
2357 //standard end of comment
2358 return 0;
2359 }
2360 */
2361
2362 TComplex FemtoFsiParsed::GetC(Int_t i) const { return TComplex(fsi_c__._2.c__[i].r, fsi_c__._2.c__[i].i); }
2363
2364 void FemtoFsiParsed::SetC(TComplex Z, Int_t i) {
2365 fsi_c__._2.c__[i].r = Z.Re();
2366 fsi_c__._2.c__[i].i = Z.Im();
2367 }
2368} // namespace Hal