62 using namespace G4InuclSpecialFunctions;
75 fissionStore.setVerboseLevel(verboseLevel);
78 getTargetData(target);
83 G4double TEM = std::sqrt(EEXS / PARA);
86 TETA = TETA / std::sinh(TETA);
107 for (
G4int i = 0; i < 50 && A1 > 30; i++) {
112 Z1 =
G4lrint(getZopt(A1, A2,
Z, X3, X4, R12) - 1.);
117 potentialMinimization(VPOT, EDEF1, VCOUL, A1, A2, Z1, Z2, AL1, BET1, R12);
129 G4double C1 = std::sqrt(getC2(A1, A2, X3, X4, R12) / TEM);
132 DZ = DZ > 0.0 ? DZ + 0.5 : -std::fabs(DZ - 0.5);
137 G4double EZ = (DMT1 + (DMT - DMT1) * TETA - VPOT + DEfin) / TEM;
139 if (EZ >= ALMA) ALMA = EZ;
140 G4double EK = VCOUL + DEfin + 0.5 * TEM;
143 if (EV > 0.0) fissionStore.addConfig(A1, Z1, EZ, EK, EV);
147 G4int store_size = fissionStore.size();
148 if (store_size == 0)
return;
151 fissionStore.generateConfiguration(ALMA,
inuclRndm());
162 G4double pmod = std::sqrt(0.001 * EK * mass1 * mass2 / mass_in);
168 G4double EV = 1000.0 * (e_in - e_out) /
A;
169 if (EV <= 0.0)
return;
185 if (verboseLevel > 3) {
189 G4double C2 = 124.57 * (1.0 / A1 + 1.0 / A2) + 0.78 * (X3 + X4) - 176.9 *
190 ((X3*X3*X3*X3) + (X4*X4*X4*X4)) + 219.36 * (1.0 / (A1 * A1) + 1.0 / (A2 * A2)) - 1.108 / R12;
202 if (verboseLevel > 3) {
206 G4double Zopt = (87.7 * (X4 - X3) * (1.0 - 1.25 * (X4 + X3)) +
207 ZT * ((124.57 / A2 + 0.78 * X4 - 176.9 * (X4*X4*X4*X4) + 219.36 / (A2 * A2)) - 0.554 / R12)) /
208 getC2(A1, A2, X3, X4, R12);
224 if (verboseLevel > 3) {
225 G4cout <<
" >>> G4Fissioner::potentialMinimization" <<
G4endl;
229 const G4int itry_max = 2000;
232 const G4double DS2 = 1.0 / DS1 / DS1;
233 G4int A1[2] = { AF, AS };
245 for (i = 0; i < 2; i++) {
248 Y2 = Z1[i] * Z1[i] / R[i];
249 C[i] = 6.8 * Y1 - 0.142 *
Y2;
250 F[i] = 12.138 * Y1 - 0.145 *
Y2;
264 while (itry < itry_max) {
268 for (i = 0; i < 2; i++) {
269 S += R[i] * (1.0 + AL1[i] + BET1[i] - 0.257 * AL1[i] * BET1[i]);
275 for (i = 0; i < 2; i++) {
276 SAL[i] = R[i] * (1.0-0.257 * BET1[i]);
277 SBE[i] = R[i] * (1.0-0.257 * AL1[i]);
280 X2[i] = X[i] * X1[i];
281 Y1 += AL1[i] * X1[i];
282 Y2 += BET1[i] * X2[i];
283 R12 += R[i] * (1.0 - AL1[i] * (1.0 - 0.6 * X[i]) + BET1[i] * (1.0 - 0.429 * X1[i]));
291 for (i = 0; i < 2; i++) {
292 RAL[i] = -R[i] * (1.0 - 0.6 * X[i]) + SAL[i] * Y3;
293 RBE[i] = R[i] * (1.0 - 0.429 * X1[i]) + SBE[i] * Y3;
299 for (i = 0; i < 2; i++) {
301 for (
G4int j = 0; j < 2; j++) {
306 if (std::fabs(AL1[i]) >= DS1) {
309 DX1 = 2.0 * (1.0 + 2.0 * AL1[i] * AL1[i] * DS2) * DEX * DS2;
312 if (std::fabs(BET1[i]) >= DS1) {
315 DX2 = 2.0 * (1.+2.0 * BET1[i] * BET1[i] * DS2) * DEX * DS2;
319 AA[i][j] = R3 * RBE[i] * RBE[j] -
322 X1[j] * SAL[i]) + SAL[i] * SAL[j] * Y4) +
323 DEL * C[i] + DEL1 * DX1;
326 AA[i1][j1] = R3 * RBE[i] * RBE[j]
329 X2[j] * SBE[i]) + SBE[i] * SBE[j] * Y4) +
330 DEL * F[i] + DEL1 * DX2;
331 AA[i][j1] = R3 * RAL[i] * RBE[j] -
334 0.6 * X1[i] * SBE[j]) + SBE[j] * SAL[i] * Y4 -
335 0.257 * R[i] * Y3 * DEL1);
336 AA[j1][i] = AA[i][j1];
340 for (i = 0; i < 2; i++) {
344 if (std::fabs(AL1[i]) >= DS1) DX1 = 2.0 * AL1[i] * DS2 *
G4Exp(AL1[i] * AL1[i] * DS2);
346 if (std::fabs(BET1[i]) >= DS1) DX2 = 2.0 * BET1[i] * DS2 *
G4Exp(BET1[i] * BET1[i] * DS2);
347 B[i] = R2 * RAL[i] - 2.0e-3 * C[i] * AL1[i] + DX1;
348 B[i + 2] = R2 * RBE[i] - 2.0e-3 * F[i] * BET1[i] + DX2;
354 for (i = 0; i < 4; i++) {
357 for (
G4int j = 0; j < 4; j++) ST1 += AA[i][j] * B[i] * B[j];
363 for (i = 0; i < 2; i++) {
364 AL1[i] += B[i] * STEP;
365 BET1[i] += B[i + 2] * STEP;
366 DSOL += B[i] * B[i] + B[i + 2] * B[i + 2];
368 DSOL = std::sqrt(DSOL);
370 if (DSOL < DSOL1)
break;
373 if (verboseLevel > 3) {
374 if (itry == itry_max)
375 G4cout <<
" maximal number of iterations in potentialMinimization " <<
G4endl
376 <<
" A1 " << AF <<
" Z1 " << ZF <<
G4endl;
380 for (i = 0; i < 2; i++) ED[i] = F[i] * BET1[i] * BET1[i] + C[i] * AL1[i] * AL1[i];
383 VP = VC + ED[0] + ED[1];