59 using namespace G4InuclSpecialFunctions;
72 fissionStore.setVerboseLevel(verboseLevel);
75 getTargetData(target);
79 G4double TEM = std::sqrt(EEXS / PARA);
82 TETA = TETA / std::sinh(TETA);
94 G4double DTEM = (A < 220 ? 0.5 : 1.15);
103 for (
G4int i = 0; i < 50 && A1 > 30; i++) {
108 Z1 =
G4lrint(getZopt(A1, A2,
Z, X3, X4, R12) - 1.);
113 potentialMinimization(VPOT, EDEF1, VCOUL, A1, A2, Z1, Z2, AL1, BET1, R12);
125 G4double C1 = std::sqrt(getC2(A1, A2, X3, X4, R12) / TEM);
128 DZ = DZ > 0.0 ? DZ + 0.5 : -std::fabs(DZ - 0.5);
133 G4double EZ = (DMT1 + (DMT - DMT1) * TETA - VPOT + DEfin) / TEM;
135 if (EZ >= ALMA) ALMA = EZ;
136 G4double EK = VCOUL + DEfin + 0.5 * TEM;
139 if (EV > 0.0) fissionStore.addConfig(A1, Z1, EZ, EK, EV);
143 G4int store_size = fissionStore.size();
144 if (store_size == 0)
return;
147 fissionStore.generateConfiguration(ALMA,
inuclRndm());
158 G4double pmod = std::sqrt(0.001 * EK * mass1 * mass2 / mass_in);
164 G4double EV = 1000.0 * (e_in - e_out) / A;
165 if (EV <= 0.0)
return;
181 if (verboseLevel > 3) {
185 G4double C2 = 124.57 * (1.0 / A1 + 1.0 / A2) + 0.78 * (X3 + X4) - 176.9 *
186 ((X3*X3*X3*X3) + (X4*X4*X4*X4)) + 219.36 * (1.0 / (A1 * A1) + 1.0 / (A2 * A2)) - 1.108 / R12;
198 if (verboseLevel > 3) {
202 G4double Zopt = (87.7 * (X4 - X3) * (1.0 - 1.25 * (X4 + X3)) +
203 ZT * ((124.57 / A2 + 0.78 * X4 - 176.9 * (X4*X4*X4*X4) + 219.36 / (A2 * A2)) - 0.554 / R12)) /
204 getC2(A1, A2, X3, X4, R12);
209 void G4Fissioner::potentialMinimization(
G4double& VP,
220 if (verboseLevel > 3) {
221 G4cout <<
" >>> G4Fissioner::potentialMinimization" <<
G4endl;
225 const G4int itry_max = 2000;
228 const G4double DS2 = 1.0 / DS1 / DS1;
229 G4int A1[2] = { AF, AS };
241 for (i = 0; i < 2; i++) {
244 Y2 = Z1[i] * Z1[i] / R[i];
245 C[i] = 6.8 * Y1 - 0.142 *
Y2;
246 F[i] = 12.138 * Y1 - 0.145 *
Y2;
260 while (itry < itry_max) {
264 for (i = 0; i < 2; i++) {
265 S += R[i] * (1.0 + AL1[i] + BET1[i] - 0.257 * AL1[i] * BET1[i]);
271 for (i = 0; i < 2; i++) {
272 SAL[i] = R[i] * (1.0-0.257 * BET1[i]);
273 SBE[i] = R[i] * (1.0-0.257 * AL1[i]);
276 X2[i] = X[i] * X1[i];
277 Y1 += AL1[i] * X1[i];
278 Y2 += BET1[i] * X2[i];
279 R12 += R[i] * (1.0 - AL1[i] * (1.0 - 0.6 * X[i]) + BET1[i] * (1.0 - 0.429 * X1[i]));
287 for (i = 0; i < 2; i++) {
288 RAL[i] = -R[i] * (1.0 - 0.6 * X[i]) + SAL[i] * Y3;
289 RBE[i] = R[i] * (1.0 - 0.429 * X1[i]) + SBE[i] * Y3;
295 for (i = 0; i < 2; i++) {
297 for (
G4int j = 0; j < 2; j++) {
302 if (std::fabs(AL1[i]) >= DS1) {
304 G4double DEX = XXX > 100.0 ? huge_num : std::exp(XXX);
305 DX1 = 2.0 * (1.0 + 2.0 * AL1[i] * AL1[i] * DS2) * DEX * DS2;
308 if (std::fabs(BET1[i]) >= DS1) {
310 G4double DEX = XXX > 100.0 ? huge_num : std::exp(XXX);
311 DX2 = 2.0 * (1.+2.0 * BET1[i] * BET1[i] * DS2) * DEX * DS2;
315 AA[i][j] = R3 * RBE[i] * RBE[j] -
318 X1[j] * SAL[i]) + SAL[i] * SAL[j] * Y4) +
319 DEL * C[i] + DEL1 * DX1;
322 AA[i1][j1] = R3 * RBE[i] * RBE[j]
325 X2[j] * SBE[i]) + SBE[i] * SBE[j] * Y4) +
326 DEL * F[i] + DEL1 * DX2;
327 AA[i][j1] = R3 * RAL[i] * RBE[j] -
330 0.6 * X1[i] * SBE[j]) + SBE[j] * SAL[i] * Y4 -
331 0.257 * R[i] * Y3 * DEL1);
332 AA[j1][i] = AA[i][j1];
336 for (i = 0; i < 2; i++) {
340 if (std::fabs(AL1[i]) >= DS1) DX1 = 2.0 * AL1[i] * DS2 * std::exp(AL1[i] * AL1[i] * DS2);
342 if (std::fabs(BET1[i]) >= DS1) DX2 = 2.0 * BET1[i] * DS2 * std::exp(BET1[i] * BET1[i] * DS2);
343 B[i] = R2 * RAL[i] - 2.0e-3 * C[i] * AL1[i] + DX1;
344 B[i + 2] = R2 * RBE[i] - 2.0e-3 * F[i] * BET1[i] + DX2;
350 for (i = 0; i < 4; i++) {
353 for (
G4int j = 0; j < 4; j++) ST1 += AA[i][j] * B[i] * B[j];
359 for (i = 0; i < 2; i++) {
360 AL1[i] += B[i] * STEP;
361 BET1[i] += B[i + 2] * STEP;
362 DSOL += B[i] * B[i] + B[i + 2] * B[i + 2];
364 DSOL = std::sqrt(DSOL);
366 if (DSOL < DSOL1)
break;
369 if (verboseLevel > 3) {
370 if (itry == itry_max)
371 G4cout <<
" maximal number of iterations in potentialMinimization " <<
G4endl
372 <<
" A1 " << AF <<
" Z1 " << ZF <<
G4endl;
376 for (i = 0; i < 2; i++) ED[i] = F[i] * BET1[i] * BET1[i] + C[i] * AL1[i] * AL1[i];
379 VP = VC + ED[0] + ED[1];
G4double randomGauss(G4double sigma)
G4double bindingEnergyAsymptotic(G4int A, G4int Z)
void setVectM(const Hep3Vector &spatial, double mass)
G4GLOB_DLL std::ostream G4cout
G4double getNucleiMass() const
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
G4double G4cbrt(G4double x)
void addRecoilFragment(const G4Fragment *aFragment)
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
G4double bindingEnergy(G4int A, G4int Z)
G4double nucleiLevelDensity(G4int A)