63 using namespace G4InuclSpecialFunctions;
76 fissionStore.setVerboseLevel(verboseLevel);
79 getTargetData(target);
84 G4double TEM = std::sqrt(EEXS / PARA);
87 TETA = TETA / std::sinh(TETA);
108 for (
G4int i = 0; i < 50 && A1 > 30; i++) {
113 Z1 =
G4lrint(getZopt(A1, A2,
Z, X3, X4, R12) - 1.);
118 potentialMinimization(VPOT, EDEF1, VCOUL, A1, A2, Z1, Z2, AL1, BET1, R12);
130 G4double C1 = std::sqrt(getC2(A1, A2, X3, X4, R12) / TEM);
133 DZ = DZ > 0.0 ? DZ + 0.5 : -std::fabs(DZ - 0.5);
138 G4double EZ = (DMT1 + (DMT - DMT1) * TETA - VPOT + DEfin) / TEM;
140 if (EZ >= ALMA) ALMA = EZ;
141 G4double EK = VCOUL + DEfin + 0.5 * TEM;
144 if (EV > 0.0) fissionStore.addConfig(A1, Z1, EZ, EK, EV);
148 G4int store_size = fissionStore.size();
149 if (store_size == 0)
return;
152 fissionStore.generateConfiguration(ALMA,
inuclRndm());
163 G4double pmod = std::sqrt(0.001 * EK * mass1 * mass2 / mass_in);
169 G4double EV = 1000.0 * (e_in - e_out) /
A;
170 if (EV <= 0.0)
return;
186 if (verboseLevel > 3) {
190 G4double C2 = 124.57 * (1.0 / A1 + 1.0 / A2) + 0.78 * (X3 + X4) - 176.9 *
191 ((X3*X3*X3*X3) + (X4*X4*X4*X4)) + 219.36 * (1.0 / (A1 * A1) + 1.0 / (A2 * A2)) - 1.108 / R12;
203 if (verboseLevel > 3) {
207 G4double Zopt = (87.7 * (X4 - X3) * (1.0 - 1.25 * (X4 + X3)) +
208 ZT * ((124.57 / A2 + 0.78 * X4 - 176.9 * (X4*X4*X4*X4) + 219.36 / (A2 * A2)) - 0.554 / R12)) /
209 getC2(A1, A2, X3, X4, R12);
225 if (verboseLevel > 3) {
226 G4cout <<
" >>> G4Fissioner::potentialMinimization" <<
G4endl;
230 const G4int itry_max = 2000;
233 const G4double DS2 = 1.0 / DS1 / DS1;
234 G4int A1[2] = { AF, AS };
246 for (i = 0; i < 2; i++) {
249 Y2 = Z1[i] * Z1[i] / R[i];
250 C[i] = 6.8 * Y1 - 0.142 *
Y2;
251 F[i] = 12.138 * Y1 - 0.145 *
Y2;
265 while (itry < itry_max) {
269 for (i = 0; i < 2; i++) {
270 S += R[i] * (1.0 + AL1[i] + BET1[i] - 0.257 * AL1[i] * BET1[i]);
276 for (i = 0; i < 2; i++) {
277 SAL[i] = R[i] * (1.0-0.257 * BET1[i]);
278 SBE[i] = R[i] * (1.0-0.257 * AL1[i]);
281 X2[i] = X[i] * X1[i];
282 Y1 += AL1[i] * X1[i];
283 Y2 += BET1[i] * X2[i];
284 R12 += R[i] * (1.0 - AL1[i] * (1.0 - 0.6 * X[i]) + BET1[i] * (1.0 - 0.429 * X1[i]));
292 for (i = 0; i < 2; i++) {
293 RAL[i] = -R[i] * (1.0 - 0.6 * X[i]) + SAL[i] * Y3;
294 RBE[i] = R[i] * (1.0 - 0.429 * X1[i]) + SBE[i] * Y3;
300 for (i = 0; i < 2; i++) {
302 for (
G4int j = 0; j < 2; j++) {
307 if (std::fabs(AL1[i]) >= DS1) {
310 DX1 = 2.0 * (1.0 + 2.0 * AL1[i] * AL1[i] * DS2) * DEX * DS2;
313 if (std::fabs(BET1[i]) >= DS1) {
316 DX2 = 2.0 * (1.+2.0 * BET1[i] * BET1[i] * DS2) * DEX * DS2;
320 AA[i][j] = R3 * RBE[i] * RBE[j] -
323 X1[j] * SAL[i]) + SAL[i] * SAL[j] * Y4) +
324 DEL * C[i] + DEL1 * DX1;
327 AA[i1][j1] = R3 * RBE[i] * RBE[j]
330 X2[j] * SBE[i]) + SBE[i] * SBE[j] * Y4) +
331 DEL * F[i] + DEL1 * DX2;
332 AA[i][j1] = R3 * RAL[i] * RBE[j] -
335 0.6 * X1[i] * SBE[j]) + SBE[j] * SAL[i] * Y4 -
336 0.257 * R[i] * Y3 * DEL1);
337 AA[j1][i] = AA[i][j1];
341 for (i = 0; i < 2; i++) {
345 if (std::fabs(AL1[i]) >= DS1) DX1 = 2.0 * AL1[i] * DS2 *
G4Exp(AL1[i] * AL1[i] * DS2);
347 if (std::fabs(BET1[i]) >= DS1) DX2 = 2.0 * BET1[i] * DS2 *
G4Exp(BET1[i] * BET1[i] * DS2);
348 B[i] = R2 * RAL[i] - 2.0e-3 * C[i] * AL1[i] + DX1;
349 B[i + 2] = R2 * RBE[i] - 2.0e-3 * F[i] * BET1[i] + DX2;
355 for (i = 0; i < 4; i++) {
358 for (
G4int j = 0; j < 4; j++) ST1 += AA[i][j] * B[i] * B[j];
364 for (i = 0; i < 2; i++) {
365 AL1[i] += B[i] * STEP;
366 BET1[i] += B[i + 2] * STEP;
367 DSOL += B[i] * B[i] + B[i + 2] * B[i + 2];
369 DSOL = std::sqrt(DSOL);
371 if (DSOL < DSOL1)
break;
374 if (verboseLevel > 3) {
375 if (itry == itry_max)
376 G4cout <<
" maximal number of iterations in potentialMinimization " <<
G4endl
377 <<
" A1 " << AF <<
" Z1 " << ZF <<
G4endl;
381 for (i = 0; i < 2; i++) ED[i] = F[i] * BET1[i] * BET1[i] + C[i] * AL1[i] * AL1[i];
384 VP = VC + ED[0] + ED[1];
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
G4double randomGauss(G4double sigma)
void addRecoilFragment(const G4Fragment *aFragment)
G4double G4cbrt(G4double x)
double A(double temperature)
G4double nucleiLevelDensity(G4int A)
G4double getZopt(G4int A1, G4int A2, G4int ZT, G4double X3, G4double X4, G4double R12) const
G4double getC2(G4int A1, G4int A2, G4double X3, G4double X4, G4double R12) const
G4double getNucleiMass() const
G4double bindingEnergyAsymptotic(G4int A, G4int Z)
G4GLOB_DLL std::ostream G4cout
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 bindingEnergy(G4int A, G4int Z)
double B(double temperature)
void setVectM(const Hep3Vector &spatial, double mass)