86 TETA = TETA / std::sinh(TETA);
107 for (
G4int i = 0; i < 50 && A1 > 30; i++) {
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;
148 if (store_size == 0)
return;
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;
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;
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);
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]));
286 G4double Y3 = -0.6 * Y1 + 0.857 * Y2;
287 G4double Y4 = (1.2 * Y1 - 2.571 * Y2) /
S;
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;
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];
G4double C(G4double temp)
G4double S(G4double temp)
G4double B(G4double temperature)
G4double D(G4double temp)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4GLOB_DLL std::ostream G4cout
void setVectM(const Hep3Vector &spatial, double mass)
void getTargetData(const G4Fragment &target)
const G4Fragment & makeFragment(G4LorentzVector mom, G4int A, G4int Z, G4double EX=0.)
void addRecoilFragment(const G4Fragment *aFragment)
void addConfig(G4double a, G4double z, G4double ez, G4double ek, G4double ev)
void setVerboseLevel(G4int verbose=1)
G4FissionConfiguration generateConfiguration(G4double amax, G4double rand) const
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
void potentialMinimization(G4double &VP, G4double(&ED)[2], G4double &VC, G4int AF, G4int AS, G4int ZF, G4int ZS, G4double AL1[2], G4double BET1[2], G4double &R12) const
G4double getZopt(G4int A1, G4int A2, G4int ZT, G4double X3, G4double X4, G4double R12) const
G4FissionStore fissionStore
G4double getC2(G4int A1, G4int A2, G4double X3, G4double X4, G4double R12) const
G4double getNucleiMass() const
G4double bindingEnergy(G4int A, G4int Z)
G4double nucleiLevelDensity(G4int A)
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
G4double randomGauss(G4double sigma)
G4double G4cbrt(G4double x)
G4double bindingEnergyAsymptotic(G4int A, G4int Z)
static const G4double Z1[5]
static const G4double AA[5]