83{
85 DefineCoefficients(pol0, pol1);
86
87
88 fDice = 0.;
90 0.125 * ((-1. /
sqr(gam + 1.)) /
sqr(eps) +
91 ((
sqr(gam) + 4. *
gam - 1.) /
sqr(gam + 1.)) / eps - 1.);
92
98
99
100 G4double helpVar2 = 0., helpVar1 = 0.;
101
102
104 helpVar2 = -1. /
sqr(gam + 1.);
106 fUnpXS = 0.125 * calcVector * sumEpsVector;
107
108
109 helpVar2 = 1. /
sqr(gam + 1.);
110 helpVar1 = -(
gam *
gam + 4. *
gam + 1.) /
sqr(gam + 1.);
111 calcVector =
G4ThreeVector(helpVar2, helpVar1, 0.5 * (gam + 3.));
112 ISPxx = 0.25 * (calcVector * sumEpsVector) / (gam - 1.);
113
114 helpVar1 = 1. /
sqr(gam + 1.);
115 calcVector =
G4ThreeVector(-helpVar1, 2. * gam * helpVar1, -1.);
116 ISPyy = 0.125 * calcVector * sumEpsVector;
117
118 helpVar1 = 1. / (
gam - 1.);
119 helpVar2 = 1. /
sqr(gam + 1.);
121 -(gam * gam + 1.) * helpVar2,
122 (gam * gam * (gam + 1.) + 7. * gam + 3.) * helpVar2, -(gam + 3.));
123 ISPzz = 0.125 * helpVar1 * (calcVector * sumEpsVector);
124
125 helpVar1 = std::sqrt(std::fabs(eps * (1. - eps) * 2. * (gam + 1.) - 1.));
126 calcVector =
G4ThreeVector(-1. / (gam * gam - 1.), 2. / (gam - 1.), 0.);
127 ISPnd = 0.125 * (calcVector * difEpsVector) * helpVar1;
128
129 fPolXS = ISPxx * polxx + ISPyy * polyy + ISPzz * polzz;
130 fPolXS += ISPnd * (polzx + polxz);
131 fPhi0 = fUnpXS + fPolXS;
132 fDice = symmXS;
133
134 if(polzz != 0.)
135 {
136 fDice *= (1. + (polzz * ISPzz / fUnpXS));
137 if(fDice < 0.)
138 fDice = 0.;
139 }
140
141 if(flag == 2)
142 {
143
144 G4double circ1 = 0., circ2 = 0., circ3 = 0.;
145 helpVar1 = 8. *
sqr(1. - eps) *
sqr(eps) * (
gam - 1.) *
sqr(gam + 1.) /
146 std::sqrt(gam * gam - 1.);
147 helpVar2 =
sqr(gam + 1.) *
sqr(eps) * (-2. * eps + 3.) -
148 (gam * gam + 3. * gam + 2.) * eps;
149 circ1 = helpVar2 +
gam;
150 circ1 /= helpVar1;
151 circ2 = helpVar2 + 1.;
152 circ2 /= helpVar1;
153 helpVar1 = std::sqrt(std::fabs(eps * (1. - eps) * 2. * (gam + 1.) - 1.));
154 helpVar1 /= std::sqrt(gam * gam - 1.);
156 circ3 = 0.125 * (calcVector * sumEpsVector) / (gam + 1.);
157 circ3 *= helpVar1;
158
159 fPhi2.setZ(circ2 * pol1.z() + circ1 * pol0.z() +
160 circ3 * (pol1.x() + pol0.x()));
161 fPhi3.setZ(-circ1 * pol1.z() - circ2 * pol0.z() -
162 circ3 * (pol1.x() + pol0.x()));
163
164
166 G4double linearZero = 0.125 * (calcVector * sumEpsVector) /
sqr(gam + 1.);
167
168
169 helpVar1 = std::sqrt(std::fabs(2. * (gam + 1.) * (1. - eps) * eps - 1.)) /
170 ((
gam + 1.) * eps * (1. - eps));
171 helpVar2 = helpVar1 * helpVar1;
172
173
174 G4double diagContrib = 0.125 * helpVar2 * (polxx + polyy - polzz);
176 0.125 * helpVar1 * (-polxz / (1. - eps) + polzx / eps);
177
178 fPhi2.setX(linearZero + diagContrib + nonDiagContrib);
179
180
181 nonDiagContrib = 0.125 * helpVar1 * (polxz / eps - polzx / (1. - eps));
182
183 fPhi3.setX(linearZero + diagContrib + nonDiagContrib);
184
185
186 helpVar1 =
187 std::sqrt(gam * gam - 1.) * (2. * (
gam + 1.) * eps * (1. - eps) - 1.);
188 helpVar1 /= 8. *
sqr(1. - eps) *
sqr(eps) *
sqr(gam + 1.) * (
gam - 1.);
189 helpVar2 = std::sqrt((gam * gam - 1.) *
190 std::fabs(2. * (gam + 1.) * eps * (1. - eps) - 1.));
191 helpVar2 /= 8. *
sqr(1. - eps) *
sqr(eps) *
sqr(gam + 1.) * (
gam - 1.);
192
193 G4double contrib21 = (-polxy + polyx) * helpVar1;
195 -(eps * (
gam + 1.) - 1.) * polyz + (eps * (gam + 1.) - gam) * polzy;
196
197 contrib32 *= helpVar2;
198 fPhi2.setY(contrib21 + contrib32);
199
200 contrib32 =
201 -(eps * (
gam + 1.) -
gam) * polyz + (eps * (gam + 1.) - 1.) * polzy;
202 contrib32 *= helpVar2;
203 fPhi3.setY(contrib21 + contrib32);
204 }
205 fPhi0 *= diffXSFactor;
206 fPhi2 *= diffXSFactor;
207 fPhi3 *= diffXSFactor;
208}