107 G4int NbLam0 = std::abs(nucleusS);
109 Ainit = -1 * nucleusA;
110 Zinit = -1 * nucleusZ;
111 Sinit = -1 * nucleusS;
115 G4int ZFP1 = 0, AFP1 = 0, AFPIMF = 0, ZFPIMF = 0, ZFP2 = 0, AFP2 = 0, SFP1 = 0, SFP2 = 0, SFPIMF = 0;
116 G4double vx_eva = 0.0, vy_eva = 0.0, vz_eva = 0.0;
117 G4double VX_PREF = 0., VY_PREF = 0., VZ_PREF = 00, VP1X, VP1Y, VP1Z, VXOUT, VYOUT, VZOUT, V_CM[3], VFP1_CM[3],
118 VFP2_CM[3], VIMF_CM[3], VX2OUT, VY2OUT, VZ2OUT;
119 G4double zf = 0.0, af = 0.0, mtota = 0.0, tkeimf = 0.0, jprf0 = 0.;
120 G4int ff = 0, afpnew = 0, zfpnew = 0, aprfp = 0, zprfp = 0, IOUNSTABLE = 0, ILOOP = 0, IEV_TAB = 0,
122 G4int fimf = 0, INMIN = 0, INMAX = 0;
124 G4int inum = eventnumber;
126 opt->optimfallowed = 1;
139 opt->nblan0 = NbLam0;
153 G4double T_init = 0., T_diff = 0., a_tilda = 0., a_tilda_BU = 0., EE_diff = 0., EINCL = 0., A_FINAL = 0.,
154 Z_FINAL = 0., E_FINAL = 0.;
156 G4double A_diff = 0., ASLOPE1, ASLOPE2, A_ACC, ABU_SLOPE, ABU_SUM = 0., AMEM = 0., ZMEM = 0., EMEM = 0., JMEM = 0.,
157 PX_BU_SUM = 0.0, PY_BU_SUM = 0.0, PZ_BU_SUM = 0.0, ETOT_SUM = 0., P_BU_SUM = 0., ZBU_SUM = 0.,
158 Z_Breakup_sum = 0., A_Breakup, Z_Breakup, N_Breakup, G_SYMM, CZ, Sigma_Z, Z_Breakup_Mean, ZTEMP = 0.,
161 G4double ETOT_PRF = 0.0, PXPRFP = 0., PYPRFP = 0., PZPRFP = 0., PPRFP = 0., VX1_BU = 0., VY1_BU = 0., VZ1_BU = 0.,
162 VBU2 = 0., GAMMA_REL = 1.0, Eexc_BU_SUM = 0., VX_BU_SUM = 0., VY_BU_SUM = 0., VZ_BU_SUM = 0.,
163 E_tot_BU = 0., EKIN_BU = 0., ZIMFBU = 0., AIMFBU = 0., ZFFBU = 0., AFFBU = 0., AFBU = 0., ZFBU = 0.,
164 EEBU = 0., TKEIMFBU = 0., vx_evabu = 0., vy_evabu = 0., vz_evabu = 0., Bvalue_BU = 0., P_BU = 0.,
165 ETOT_BU = 1., PX_BU = 0., PY_BU = 0., PZ_BU = 0., VX2_BU = 0., VY2_BU = 0., VZ2_BU = 0.;
167 G4int ABU_DIFF, ZBU_DIFF, NBU_DIFF;
168 G4int INEWLOOP = 0, ILOOPBU = 0;
176 std::cout <<
"Error - Remnant with a mass number A below 1." << std::endl;
181 for (
G4int j = 0; j < 3; j++)
191 for (
G4int I2 = 0; I2 < 12; I2++)
192 BU_TAB[I1][I2] = 0.0;
193 for (
G4int I2 = 0; I2 < 6; I2++)
195 BU_TAB_TEMP[I1][I2] = 0.0;
196 BU_TAB_TEMP1[I1][I2] = 0.0;
197 EV_TAB_TEMP[I1][I2] = 0.0;
198 EV_TAB[I1][I2] = 0.0;
199 EV_TAB_SSC[I1][I2] = 0.0;
200 EV_TEMP[I1][I2] = 0.0;
228 G4double pincl = std::sqrt(pxrem * pxrem + pyrem * pyrem + pzrem * pzrem);
230 G4double ETOT_incl = std::sqrt(pincl * pincl + (AAINCL * amu) * (AAINCL * amu));
231 G4double VX_incl = C * pxrem / ETOT_incl;
232 G4double VY_incl = C * pyrem / ETOT_incl;
233 G4double VZ_incl = C * pzrem / ETOT_incl;
263 if (T_freeze_out_in >= 0.0)
265 T_freeze_out = T_freeze_out_in;
269 T_freeze_out =
max(9.33 * std::exp(-0.00282 * AAINCL), 5.5);
275 a_tilda = ald->av * aprf + ald->as * std::pow(aprf, 2.0 / 3.0) + ald->ak * std::pow(aprf, 1.0 / 3.0);
277 T_init = std::sqrt(EINCL / a_tilda);
279 T_diff = T_init - T_freeze_out;
281 if (T_diff > 0.1 && zprf > 2. && (aprf - zprf) > 0.)
287 for (
G4int i = 0; i < 5; i++)
289 EE_diff = EINCL - a_tilda * T_freeze_out * T_freeze_out;
296 A_diff =
dint(EE_diff / (8.0 * 5.0 / T_freeze_out));
301 A_FINAL = AAINCL - A_diff;
304 ald->av * A_FINAL + ald->as * std::pow(A_FINAL, 2.0 / 3.0) + ald->ak * std::pow(A_FINAL, 1.0 / 3.0);
305 E_FINAL = a_tilda * T_freeze_out * T_freeze_out;
309 EE_diff = EINCL - E_FINAL;
321 Z_FINAL =
dint(zprf * A_FINAL / (aprf));
330 A_diff = AAINCL - aprf;
341 else if (A_diff > 1.0)
350 a_tilda = ald->av * AAINCL + ald->as * std::pow(AAINCL, 2.0 / 3.0) + ald->ak * std::pow(AAINCL, 1.0 / 3.0);
352 E_FINAL = a_tilda * T_freeze_out * T_freeze_out;
354 ABU_SLOPE = (ASLOPE1 - ASLOPE2) / 4.0 * (E_FINAL / AAINCL) + ASLOPE1 - (ASLOPE1 - ASLOPE2) * 7.0 / 4.0;
366 if (ABU_SLOPE > -1.01)
370 Z_Breakup_sum = Z_FINAL;
374 for (
G4int i = 0; i < 100; i++)
383 std::cout <<
"WARNING: IPOWERLIMHAZ CALLED MORE THAN 100 TIMES WHEN "
384 "CALCULATING A_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED: "
385 << A_Breakup << std::endl;
389 if (A_Breakup > AAINCL)
392 if (A_Breakup <= 0.0)
394 std::cout <<
"A_BREAKUP <= 0 " << std::endl;
398 A_ACC = A_ACC + A_Breakup;
403 Z_Breakup_Mean =
dint(A_Breakup * ZAINCL / AAINCL);
405 Z_Breakup_sum = Z_Breakup_sum + Z_Breakup_Mean;
408 G_SYMM = 34.2281 - 5.14037 * E_FINAL / AAINCL;
409 if (E_FINAL / AAINCL < 2.0)
411 if (E_FINAL / AAINCL > 4.0)
417 CZ = 2.0 * G_SYMM * 4.0 / A_Breakup;
421 Sigma_Z = std::sqrt(T_freeze_out / CZ);
430 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN "
431 "CALCULATING Z_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE "
433 << A_Breakup <<
" " << Z_Breakup << std::endl;
439 if ((A_Breakup - Z_Breakup) < 0.0)
441 if ((A_Breakup - Z_Breakup) == 0.0 && Z_Breakup != 1.0)
444 if (Z_Breakup >= ZAINCL)
449 std::cout <<
"Z_BREAKUP RESAMPLED MORE THAN 10 TIMES; EVENT WILL "
450 "BE RESAMPLED AGAIN "
462 if (
idnint(A_Breakup - Z_Breakup) < INMIN ||
idnint(A_Breakup - Z_Breakup) > (INMAX + 5))
475 N_Breakup = A_Breakup - Z_Breakup;
476 BU_TAB[I_Breakup][0] =
dint(Z_Breakup);
477 BU_TAB[I_Breakup][1] =
dint(A_Breakup);
478 ABU_SUM = ABU_SUM + BU_TAB[i][1];
479 ZBU_SUM = ZBU_SUM + BU_TAB[i][0];
482 BU_TAB[I_Breakup][3] = 0.0;
483 I_Breakup = I_Breakup + 1;
484 IMULTBU = IMULTBU + 1;
504 ABU_DIFF =
idnint(ABU_SUM + aprf - AAINCL);
505 ZBU_DIFF =
idnint(ZBU_SUM + zprf - ZAINCL);
506 NBU_DIFF =
idnint((ABU_SUM - ZBU_SUM) + (aprf - zprf) - (AAINCL - ZAINCL));
509 std::cout <<
"WARNING - MORE THAN 200 BU " << IMULTBU << std::endl;
512 std::cout <<
"WARNING - LESS THAN 1 BU " << IMULTBU << std::endl;
516 for (
G4int i = 0; i < IMULTBU; i++)
519 while (NBU_DIFF != 0 && ZBU_DIFF != 0)
531 std::cout <<
"WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING "
532 "N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED."
537 if (IMEM_BU[IEL] == 1)
540 std::cout <<
"5555:" << IEL << RHAZ << IMULTBU << std::endl;
542 std::cout <<
"5555:" << IEL << RHAZ << IMULTBU << std::endl;
545 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0] -
DSIGN(1.0,
G4double(NBU_DIFF)));
547 else if (IEL > IMULTBU)
560 else if (IEL > IMULTBU)
569 if (ZTEMP < 1.0 && N_Breakup < 1.0)
583 BU_TAB[IEL][0] =
dint(ZTEMP);
584 BU_TAB[IEL][1] =
dint(ZTEMP + N_Breakup);
586 else if (IEL > IMULTBU)
589 aprf =
dint(ZTEMP + N_Breakup);
591 NBU_DIFF = NBU_DIFF -
ISIGN(1, NBU_DIFF);
592 ZBU_DIFF = ZBU_DIFF -
ISIGN(1, ZBU_DIFF);
596 for (
G4int i = 0; i < IMULTBU; i++)
599 if (NBU_DIFF != 0 && ZBU_DIFF == 0)
601 while (NBU_DIFF > 0 || NBU_DIFF < 0)
609 std::cout <<
"WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING "
610 "N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED."
615 if (IMEM_BU[IEL] == 1)
618 if (IPROBA > IMULTBU + 1 && NBU_DIFF > 0)
620 std::cout <<
"###',IPROBA,IMULTBU,NBU_DIFF,ZBU_DIFF,T_freeze_out" << std::endl;
624 BU_TAB[IEL][1] =
dint(BU_TAB[IEL][1] -
G4double(NBU_DIFF));
634 std::cout <<
"5556:" << IEL << RHAZ << IMULTBU << std::endl;
636 std::cout <<
"5556:" << IEL << RHAZ << IMULTBU << std::endl;
639 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0] -
DSIGN(1.0,
G4double(NBU_DIFF)));
641 else if (IEL > IMULTBU)
652 ATEMP =
dint(BU_TAB[IEL][0] + N_Breakup);
654 else if (IEL > IMULTBU)
656 ATEMP =
dint(zprf + N_Breakup);
658 if ((ATEMP - N_Breakup) < 1.0 && N_Breakup < 1.0)
669 BU_TAB[IEL][1] =
dint(BU_TAB[IEL][0] + N_Breakup);
670 else if (IEL > IMULTBU)
671 aprf =
dint(zprf + N_Breakup);
673 NBU_DIFF = NBU_DIFF -
ISIGN(1, NBU_DIFF);
677 for (
G4int i = 0; i < IMULTBU; i++)
682 if (ZBU_DIFF != 0 && NBU_DIFF == 0)
684 while (ZBU_DIFF > 0 || ZBU_DIFF < 0)
692 std::cout <<
"WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING "
693 "N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED."
698 if (IMEM_BU[IEL] == 1)
701 if (IPROBA > IMULTBU + 1 && ZBU_DIFF > 0)
703 std::cout <<
"###',IPROBA,IMULTBU,NBU_DIFF,ZBU_DIFF,T_freeze_out" << std::endl;
707 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0]);
708 BU_TAB[IEL][0] =
dint(BU_TAB[IEL][0] -
G4double(ZBU_DIFF));
709 BU_TAB[IEL][1] =
dint(BU_TAB[IEL][0] + N_Breakup);
715 N_Breakup = aprf - zprf;
717 aprf =
dint(zprf + N_Breakup);
723 std::cout <<
"5557:" << IEL << RHAZ << IMULTBU << std::endl;
725 std::cout <<
"5557:" << IEL << RHAZ << IMULTBU << std::endl;
728 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0]);
731 else if (IEL > IMULTBU)
733 N_Breakup =
dint(aprf - zprf);
736 ATEMP =
dint(ZTEMP + N_Breakup);
742 if ((ATEMP - ZTEMP) < 0.0)
747 if ((ATEMP - ZTEMP) < 1.0 && ZTEMP < 1.0)
754 BU_TAB[IEL][0] =
dint(ZTEMP);
755 BU_TAB[IEL][1] =
dint(ZTEMP + N_Breakup);
762 aprf =
dint(ZTEMP + N_Breakup);
765 ZBU_DIFF = ZBU_DIFF -
ISIGN(1, ZBU_DIFF);
775 for (
G4int i = 0; i < IMULTBU; i++)
782 if (BU_TAB[i][0] > 2.0)
784 a_tilda_BU = ald->av * BU_TAB[i][1] + ald->as * std::pow(BU_TAB[i][1], 2.0 / 3.0) +
785 ald->ak * std::pow(BU_TAB[i][1], 1.0 / 3.0);
786 BU_TAB[i][2] = a_tilda_BU * T_freeze_out * T_freeze_out;
793 if (BU_TAB[i][0] > ZMEM)
805 BU_TAB[IMEM][0] = zprf;
806 BU_TAB[IMEM][1] = aprf;
807 BU_TAB[IMEM][2] = ee;
808 BU_TAB[IMEM][3] = jprf;
820 for (
G4int i = 0; i < IMULTBU; i++)
822 ABU_SUM = ABU_SUM + BU_TAB[i][1];
823 ZBU_SUM = ZBU_SUM + BU_TAB[i][0];
825 ABU_DIFF =
idnint(ABU_SUM - AAINCL);
826 ZBU_DIFF =
idnint(ZBU_SUM - ZAINCL);
828 if (ABU_DIFF != 0 || ZBU_DIFF != 0)
829 std::cout <<
"Problem of mass in BU " << ABU_DIFF <<
" " << ZBU_DIFF << std::endl;
836 AMOMENT(AAINCL, aprf, 1, &PXPRFP, &PYPRFP, &PZPRFP);
837 PPRFP = std::sqrt(PXPRFP * PXPRFP + PYPRFP * PYPRFP + PZPRFP * PZPRFP);
840 ETOT_PRF = std::sqrt(PPRFP * PPRFP + (aprf * amu) * (aprf * amu));
841 VX_PREF = C * PXPRFP / ETOT_PRF;
842 VY_PREF = C * PYPRFP / ETOT_PRF;
843 VZ_PREF = C * PZPRFP / ETOT_PRF;
846 tke_bu(zprf, aprf, ZAINCL, AAINCL, &VX1_BU, &VY1_BU, &VZ1_BU);
853 lorentz_boost(VX1_BU, VY1_BU, VZ1_BU, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT, &VZOUT);
860 VBU2 = VX_PREF * VX_PREF + VY_PREF * VY_PREF + VZ_PREF * VZ_PREF;
861 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C * C));
862 ETOT_PRF = aprf * amu / GAMMA_REL;
863 PXPRFP = ETOT_PRF * VX_PREF / C;
864 PYPRFP = ETOT_PRF * VY_PREF / C;
865 PZPRFP = ETOT_PRF * VZ_PREF / C;
879 for (I_Breakup = 0; I_Breakup < IMULTBU; I_Breakup++)
882 Bvalue_BU = Bvalue_BU +
eflmac(
idnint(BU_TAB[I_Breakup][1]),
idnint(BU_TAB[I_Breakup][0]), 1, 0);
883 Eexc_BU_SUM = Eexc_BU_SUM + BU_TAB[I_Breakup][2];
885 AMOMENT(AAINCL, BU_TAB[I_Breakup][1], 1, &PX_BU, &PY_BU, &PZ_BU);
886 P_BU = std::sqrt(PX_BU * PX_BU + PY_BU * PY_BU + PZ_BU * PZ_BU);
889 ETOT_BU = std::sqrt(P_BU * P_BU + (BU_TAB[I_Breakup][1] * amu) * (BU_TAB[I_Breakup][1] * amu));
890 BU_TAB[I_Breakup][4] = C * PX_BU / ETOT_BU;
891 BU_TAB[I_Breakup][5] = C * PY_BU / ETOT_BU;
892 BU_TAB[I_Breakup][6] = C * PZ_BU / ETOT_BU;
894 tke_bu(BU_TAB[I_Breakup][0], BU_TAB[I_Breakup][1], ZAINCL, AAINCL, &VX2_BU, &VY2_BU, &VZ2_BU);
904 BU_TAB[I_Breakup][4],
905 BU_TAB[I_Breakup][5],
906 BU_TAB[I_Breakup][6],
911 BU_TAB[I_Breakup][4] = VXOUT;
912 BU_TAB[I_Breakup][5] = VYOUT;
913 BU_TAB[I_Breakup][6] = VZOUT;
916 VBU2 = BU_TAB[I_Breakup][4] * BU_TAB[I_Breakup][4] + BU_TAB[I_Breakup][5] * BU_TAB[I_Breakup][5] +
917 BU_TAB[I_Breakup][6] * BU_TAB[I_Breakup][6];
918 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C * C));
919 ETOT_BU = BU_TAB[I_Breakup][1] * amu / GAMMA_REL;
920 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] / C;
921 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] / C;
922 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] / C;
924 PX_BU_SUM = PX_BU_SUM + PX_BU;
925 PY_BU_SUM = PY_BU_SUM + PY_BU;
926 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
931 P_BU_SUM = std::sqrt(PX_BU_SUM * PX_BU_SUM + PY_BU_SUM * PY_BU_SUM + PZ_BU_SUM * PZ_BU_SUM);
934 ETOT_SUM = std::sqrt(P_BU_SUM * P_BU_SUM + (AAINCL * amu) * (AAINCL * amu));
936 VX_BU_SUM = C * PX_BU_SUM / ETOT_SUM;
937 VY_BU_SUM = C * PY_BU_SUM / ETOT_SUM;
938 VZ_BU_SUM = C * PZ_BU_SUM / ETOT_SUM;
945 lorentz_boost(-VX_BU_SUM, -VY_BU_SUM, -VZ_BU_SUM, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT, &VZOUT);
951 VBU2 = VX_PREF * VX_PREF + VY_PREF * VY_PREF + VZ_PREF * VZ_PREF;
952 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C * C));
953 ETOT_PRF = aprf * amu / GAMMA_REL;
954 PXPRFP = ETOT_PRF * VX_PREF / C;
955 PYPRFP = ETOT_PRF * VY_PREF / C;
956 PZPRFP = ETOT_PRF * VZ_PREF / C;
967 EKIN_BU = aprf * amu / GAMMA_REL - aprf * amu;
969 for (I_Breakup = 0; I_Breakup < IMULTBU; I_Breakup++)
979 BU_TAB[I_Breakup][4],
980 BU_TAB[I_Breakup][5],
981 BU_TAB[I_Breakup][6],
986 BU_TAB[I_Breakup][4] = VXOUT;
987 BU_TAB[I_Breakup][5] = VYOUT;
988 BU_TAB[I_Breakup][6] = VZOUT;
990 VBU2 = BU_TAB[I_Breakup][4] * BU_TAB[I_Breakup][4] + BU_TAB[I_Breakup][5] * BU_TAB[I_Breakup][5] +
991 BU_TAB[I_Breakup][6] * BU_TAB[I_Breakup][6];
992 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C * C));
994 ETOT_BU = BU_TAB[I_Breakup][1] * amu / GAMMA_REL;
996 EKIN_BU = EKIN_BU + BU_TAB[I_Breakup][1] * amu / GAMMA_REL - BU_TAB[I_Breakup][1] * amu;
998 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] / C;
999 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] / C;
1000 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] / C;
1001 E_tot_BU = E_tot_BU + ETOT_BU;
1003 PX_BU_SUM = PX_BU_SUM + PX_BU;
1004 PY_BU_SUM = PY_BU_SUM + PY_BU;
1005 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
1008 if (std::abs(PX_BU_SUM) > 10. || std::abs(PY_BU_SUM) > 10. || std::abs(PZ_BU_SUM) > 10.)
1012 P_BU_SUM = std::sqrt(PX_BU_SUM * PX_BU_SUM + PY_BU_SUM * PY_BU_SUM + PZ_BU_SUM * PZ_BU_SUM);
1015 ETOT_SUM = std::sqrt(P_BU_SUM * P_BU_SUM + (AAINCL * amu) * (AAINCL * amu));
1017 VX_BU_SUM = C * PX_BU_SUM / ETOT_SUM;
1018 VY_BU_SUM = C * PY_BU_SUM / ETOT_SUM;
1019 VZ_BU_SUM = C * PZ_BU_SUM / ETOT_SUM;
1026 lorentz_boost(-VX_BU_SUM, -VY_BU_SUM, -VZ_BU_SUM, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT, &VZOUT);
1032 VBU2 = VX_PREF * VX_PREF + VY_PREF * VY_PREF + VZ_PREF * VZ_PREF;
1033 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C * C));
1034 ETOT_PRF = aprf * amu / GAMMA_REL;
1035 PXPRFP = ETOT_PRF * VX_PREF / C;
1036 PYPRFP = ETOT_PRF * VY_PREF / C;
1037 PZPRFP = ETOT_PRF * VZ_PREF / C;
1046 E_tot_BU = ETOT_PRF;
1048 EKIN_BU = aprf * amu / GAMMA_REL - aprf * amu;
1050 for (I_Breakup = 0; I_Breakup < IMULTBU; I_Breakup++)
1060 BU_TAB[I_Breakup][4],
1061 BU_TAB[I_Breakup][5],
1062 BU_TAB[I_Breakup][6],
1067 BU_TAB[I_Breakup][4] = VXOUT;
1068 BU_TAB[I_Breakup][5] = VYOUT;
1069 BU_TAB[I_Breakup][6] = VZOUT;
1071 VBU2 = BU_TAB[I_Breakup][4] * BU_TAB[I_Breakup][4] + BU_TAB[I_Breakup][5] * BU_TAB[I_Breakup][5] +
1072 BU_TAB[I_Breakup][6] * BU_TAB[I_Breakup][6];
1073 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C * C));
1075 ETOT_BU = BU_TAB[I_Breakup][1] * amu / GAMMA_REL;
1077 EKIN_BU = EKIN_BU + BU_TAB[I_Breakup][1] * amu / GAMMA_REL - BU_TAB[I_Breakup][1] * amu;
1079 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] / C;
1080 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] / C;
1081 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] / C;
1082 E_tot_BU = E_tot_BU + ETOT_BU;
1084 PX_BU_SUM = PX_BU_SUM + PX_BU;
1085 PY_BU_SUM = PY_BU_SUM + PY_BU;
1086 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
1095 for (
G4int i = 0; i < IMULTBU; i++)
1097 if (BU_TAB[i][0] < 3.0 || BU_TAB[i][0] == BU_TAB[i][1])
1118 BU_TAB[i][4] = VP1X;
1119 BU_TAB[i][5] = VP1Y;
1120 BU_TAB[i][6] = VP1Z;
1123 for (
int IJ = 0; IJ < ILOOP; IJ++)
1125 BU_TAB[IMULTBU + INEWLOOP + IJ][0] = BU_TAB_TEMP[IJ][0];
1126 BU_TAB[IMULTBU + INEWLOOP + IJ][1] = BU_TAB_TEMP[IJ][1];
1127 BU_TAB[IMULTBU + INEWLOOP + IJ][4] = BU_TAB_TEMP[IJ][2];
1128 BU_TAB[IMULTBU + INEWLOOP + IJ][5] = BU_TAB_TEMP[IJ][3];
1129 BU_TAB[IMULTBU + INEWLOOP + IJ][6] = BU_TAB_TEMP[IJ][4];
1130 BU_TAB[IMULTBU + INEWLOOP + IJ][2] = 0.0;
1131 BU_TAB[IMULTBU + INEWLOOP + IJ][3] = 0.0;
1134 INEWLOOP = INEWLOOP + ILOOP;
1141 IMULTBU = IMULTBU + INEWLOOP;
1143 opt->optimfallowed = 1;
1152 for (
G4int i = 0; i < IMULTBU; i++)
1153 sumN = sumN + BU_TAB[i][1] - BU_TAB[i][0];
1155 for (
G4int i = 0; i < IMULTBU; i++)
1157 problamb[i] = (BU_TAB[i][1] - BU_TAB[i][0]) / sumN;
1160 Nblamb =
new G4int[IMULTBU];
1161 for (
G4int i = 0; i < IMULTBU; i++)
1163 for (
G4int j = 0; j < NbLam0;)
1165 G4double probtotal = (aprf - zprf) / sumN;
1168 if (ran <= probtotal)
1173 for (
G4int i = 0; i < IMULTBU; i++)
1176 if (probtotal < ran && ran <= probtotal + problamb[i])
1178 Nblamb[i] = Nblamb[i] + 1;
1181 probtotal = probtotal + problamb[i];
1187 for (
G4int i = 0; i < IMULTBU; i++)
1189 EEBU = BU_TAB[i][2];
1190 BU_TAB[i][10] = BU_TAB[i][6];
1192 if (BU_TAB[i][0] > 2.0)
1194 G4int nbl = Nblamb[i];
1218 BU_TAB[i][9] = jprfbu;
1222 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
1224 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1225 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1226 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1241 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1242 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1243 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1245 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1254 vx_evabu, vy_evabu, vz_evabu, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VXOUT, &VYOUT, &VZOUT);
1255 BU_TAB[i][4] = VXOUT;
1256 BU_TAB[i][5] = VYOUT;
1257 BU_TAB[i][6] = VZOUT;
1261 BU_TAB[i][7] =
dint(ZFBU);
1262 BU_TAB[i][8] =
dint(AFBU);
1263 BU_TAB[i][11] = nbl;
1274 opt->optimfallowed = 0;
1277 G4double EkinR1 = TKEIMFBU * AIMFBU / (AFBU + AIMFBU);
1278 G4double EkinR2 = TKEIMFBU * AFBU / (AFBU + AIMFBU);
1279 G4double V1 = std::sqrt(EkinR1 / AFBU) * 1.3887;
1280 G4double V2 = std::sqrt(EkinR2 / AIMFBU) * 1.3887;
1282 G4double VPERP1 = std::sqrt(
V1 *
V1 - VZ1_IMF * VZ1_IMF);
1284 G4double VX1_IMF = VPERP1 * std::sin(ALPHA1);
1285 G4double VY1_IMF = VPERP1 * std::cos(ALPHA1);
1290 G4double EEIMFP = EEBU * AFBU / (AFBU + AIMFBU);
1291 G4double EEIMF = EEBU * AIMFBU / (AFBU + AIMFBU);
1295 0.40 * 931.490 * 1.160 * 1.160 * (std::pow(AIMFBU, 5.0 / 3.0) + std::pow(AFBU, 5.0 / 3.0)) +
1296 931.490 * 1.160 * 1.160 * AIMFBU * AFBU / (AIMFBU + AFBU) *
1297 (std::pow(AIMFBU, 1. / 3.) + std::pow(AFBU, 1. / 3.)) *
1298 (std::pow(AIMFBU, 1. / 3.) + std::pow(AFBU, 1. / 3.));
1301 BU_TAB[i][9] * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(AFBU, 5.0 / 3.0) / IINERTTOT;
1303 BU_TAB[i][9] * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(AIMFBU, 5.0 / 3.0) / IINERTTOT;
1311 VX1_IMF, VY1_IMF, VZ1_IMF, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VXOUT, &VYOUT, &VZOUT);
1312 BU_TAB[i][4] = VXOUT;
1313 BU_TAB[i][5] = VYOUT;
1314 BU_TAB[i][6] = VZOUT;
1316 G4double vx1ev_imf = 0., vy1ev_imf = 0., vz1ev_imf = 0., zdummy = 0., adummy = 0., tkedummy = 0.,
1322 G4double pbH = (AFBU - ZFBU) / (AFBU - ZFBU + AIMFBU - ZIMFBU);
1323 for (
G4int j = 0; j < nbl; j++)
1357 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
1359 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1360 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1361 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1376 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1377 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1378 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1380 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1382 BU_TAB[i][7] =
dint(ZFFBU);
1383 BU_TAB[i][8] =
dint(AFFBU);
1384 BU_TAB[i][11] = NbLamH;
1398 BU_TAB[i][4] = VXOUT;
1399 BU_TAB[i][5] = VYOUT;
1400 BU_TAB[i][6] = VZOUT;
1404 opt->optimfallowed = 0;
1407 G4double zffimf, affimf, zdummy1, adummy1, tkedummy1, jprf2, vx2ev_imf, vy2ev_imf, vz2ev_imf;
1431 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
1433 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1434 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1435 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1451 lorentz_boost(VX2_IMF, VY2_IMF, VZ2_IMF, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1452 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
1453 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
1454 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
1456 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1458 BU_TAB[IMULTBU + ILOOPBU][0] = BU_TAB[i][0];
1459 BU_TAB[IMULTBU + ILOOPBU][1] = BU_TAB[i][1];
1460 BU_TAB[IMULTBU + ILOOPBU][2] = BU_TAB[i][2];
1461 BU_TAB[IMULTBU + ILOOPBU][3] = BU_TAB[i][3];
1462 BU_TAB[IMULTBU + ILOOPBU][7] =
dint(zffimf);
1463 BU_TAB[IMULTBU + ILOOPBU][8] =
dint(affimf);
1464 BU_TAB[IMULTBU + ILOOPBU][11] = NbLamimf;
1467 VX2_IMF, VY2_IMF, VZ2_IMF, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VXOUT, &VYOUT, &VZOUT);
1468 lorentz_boost(vx2ev_imf, vy2ev_imf, vz2ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1469 BU_TAB[IMULTBU + ILOOPBU][4] = VX2OUT;
1470 BU_TAB[IMULTBU + ILOOPBU][5] = VY2OUT;
1471 BU_TAB[IMULTBU + ILOOPBU][6] = VZ2OUT;
1472 ILOOPBU = ILOOPBU + 1;
1481 BU_TAB[i][7] = BU_TAB[i][0];
1482 BU_TAB[i][8] = BU_TAB[i][1];
1486 BU_TAB[i][11] = Nblamb[i];
1490 IMULTBU = IMULTBU + ILOOPBU;
1498 for (
G4int i = 0; i < IMULTBU; i++)
1500 ABU_SUM = ABU_SUM + BU_TAB[i][8];
1501 ZBU_SUM = ZBU_SUM + BU_TAB[i][7];
1524 ABU_SUM = ABU_SUM +
G4double(afpnew) - BU_TAB[i][8];
1525 ZBU_SUM = ZBU_SUM +
G4double(zfpnew) - BU_TAB[i][7];
1528 BU_TAB[i][4] = VP1X;
1529 BU_TAB[i][5] = VP1Y;
1530 BU_TAB[i][6] = VP1Z;
1533 for (
G4int IJ = 0; IJ < ILOOP; IJ++)
1535 BU_TAB[IMULTBU + INEWLOOP + IJ][7] = BU_TAB_TEMP1[IJ][0];
1536 BU_TAB[IMULTBU + INEWLOOP + IJ][8] = BU_TAB_TEMP1[IJ][1];
1537 BU_TAB[IMULTBU + INEWLOOP + IJ][4] = BU_TAB_TEMP1[IJ][2];
1538 BU_TAB[IMULTBU + INEWLOOP + IJ][5] = BU_TAB_TEMP1[IJ][3];
1539 BU_TAB[IMULTBU + INEWLOOP + IJ][6] = BU_TAB_TEMP1[IJ][4];
1540 BU_TAB[IMULTBU + INEWLOOP + IJ][2] = 0.0;
1541 BU_TAB[IMULTBU + INEWLOOP + IJ][3] = 0.0;
1542 BU_TAB[IMULTBU + INEWLOOP + IJ][0] = BU_TAB[i][0];
1543 BU_TAB[IMULTBU + INEWLOOP + IJ][1] = BU_TAB[i][1];
1544 BU_TAB[IMULTBU + INEWLOOP + IJ][11] = BU_TAB[i][11];
1545 ABU_SUM = ABU_SUM + BU_TAB[IMULTBU + INEWLOOP + IJ][8];
1546 ZBU_SUM = ZBU_SUM + BU_TAB[IMULTBU + INEWLOOP + IJ][7];
1549 INEWLOOP = INEWLOOP + ILOOP;
1554 IMULTBU = IMULTBU + INEWLOOP;
1557 lorentz_boost(VX_incl, VY_incl, VZ_incl, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT, &VZOUT);
1562 for (
G4int i = 0; i < IMULTBU; i++)
1564 lorentz_boost(VX_incl, VY_incl, VZ_incl, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VXOUT, &VYOUT, &VZOUT);
1565 BU_TAB[i][4] = VXOUT;
1566 BU_TAB[i][5] = VYOUT;
1567 BU_TAB[i][6] = VZOUT;
1569 for (
G4int i = 0; i < IEV_TAB; i++)
1571 lorentz_boost(VX_incl, VY_incl, VZ_incl, EV_TAB[i][2], EV_TAB[i][3], EV_TAB[i][4], &VXOUT, &VYOUT, &VZOUT);
1572 EV_TAB[i][2] = VXOUT;
1573 EV_TAB[i][3] = VYOUT;
1574 EV_TAB[i][4] = VZOUT;
1577 std::cout <<
"IMULTBU>200 " << IMULTBU << std::endl;
1604 opt->optimfallowed = 1;
1611 if (zprfp <= 2 && zprfp < aprfp)
1630 if (zprfp <= 2 && zprfp == aprfp)
1650 for (
G4int I = 0; I < ILOOP; I++)
1652 for (
G4int IJ = 0; IJ < 6; IJ++)
1653 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1655 IEV_TAB = IEV_TAB + ILOOP;
1691 for (
G4int I = 0; I < ILOOP; I++)
1693 for (
G4int IJ = 0; IJ < 6; IJ++)
1694 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1696 IEV_TAB = IEV_TAB + ILOOP;
1733 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
1735 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1736 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1737 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1744 VX_PREF, VY_PREF, VZ_PREF, EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1745 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1746 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1747 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1749 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1754 lorentz_boost(VX_PREF, VY_PREF, VZ_PREF, vx_eva, vy_eva, vz_eva, &VXOUT, &VYOUT, &VZOUT);
1759 if (ff == 0 && fimf == 0)
1772 VFP1_CM[0] = V_CM[0];
1773 VFP1_CM[1] = V_CM[1];
1774 VFP1_CM[2] = V_CM[2];
1775 for (
G4int j = 0; j < 3; j++)
1782 if (ff == 1 && fimf == 0)
1784 if (ff == 0 && fimf == 1)
1800 G4int IEV_TAB_FIS = 0, imode = 0;
1802 G4double vx1_fission = 0., vy1_fission = 0., vz1_fission = 0.;
1803 G4double vx2_fission = 0., vy2_fission = 0., vz2_fission = 0.;
1804 G4double vx_eva_sc = 0., vy_eva_sc = 0., vz_eva_sc = 0.;
1830 for (
G4int IJ = 0; IJ < IEV_TAB_FIS; IJ++)
1832 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1833 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1834 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1841 V_CM[0], V_CM[1], V_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1842 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1843 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1844 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1846 IEV_TAB = IEV_TAB + IEV_TAB_FIS;
1862 lorentz_boost(vx1_fission, vy1_fission, vz1_fission, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1863 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1864 VFP1_CM[0] = VX2OUT;
1865 VFP1_CM[1] = VY2OUT;
1866 VFP1_CM[2] = VZ2OUT;
1873 lorentz_boost(vx2_fission, vy2_fission, vz2_fission, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1874 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1875 VFP2_CM[0] = VX2OUT;
1876 VFP2_CM[1] = VY2OUT;
1877 VFP2_CM[2] = VZ2OUT;
1883 else if (ftype == 2)
1889 opt->optimfallowed = 1;
1894 G4double pbH = (af - zf) / (af - zf + aimf - zimf);
1896 for (
G4int i = 0; i < NbLam0; i++)
1909 G4double EkinR1 = tkeimf * aimf / (af + aimf);
1910 G4double EkinR2 = tkeimf * af / (af + aimf);
1911 G4double V1 = std::sqrt(EkinR1 / af) * 1.3887;
1912 G4double V2 = std::sqrt(EkinR2 / aimf) * 1.3887;
1914 G4double VPERP1 = std::sqrt(
V1 *
V1 - VZ1_IMF * VZ1_IMF);
1916 G4double VX1_IMF = VPERP1 * std::sin(ALPHA1);
1917 G4double VY1_IMF = VPERP1 * std::cos(ALPHA1);
1922 G4double EEIMFP = ee * af / (af + aimf);
1923 G4double EEIMF = ee * aimf / (af + aimf);
1926 G4double IINERTTOT = 0.40 * 931.490 * 1.160 * 1.160 * (std::pow(aimf, 5.0 / 3.0) + std::pow(af, 5.0 / 3.0)) +
1927 931.490 * 1.160 * 1.160 * aimf * af / (aimf + af) *
1928 (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.)) *
1929 (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.));
1931 G4double JPRFHEAVY = jprf0 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(af, 5.0 / 3.0) / IINERTTOT;
1932 G4double JPRFLIGHT = jprf0 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(aimf, 5.0 / 3.0) / IINERTTOT;
1934 std::cout <<
"RN117-4,AF,ZF,EE,JPRFheavy" << std::endl;
1936 G4double vx1ev_imf = 0., vy1ev_imf = 0., vz1ev_imf = 0., zdummy = 0., adummy = 0., tkedummy = 0., jprf1 = 0.;
1960 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
1962 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1963 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1964 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1971 V_CM[0], V_CM[1], V_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1972 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1973 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
1974 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
1975 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
1977 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1982 opt->optimfallowed = 0;
1986 G4double zffimf, affimf, zdummy1 = 0., adummy1 = 0., tkedummy1 = 0., jprf2, vx2ev_imf, vy2ev_imf, vz2ev_imf;
2010 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
2012 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
2013 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
2014 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
2021 V_CM[0], V_CM[1], V_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
2022 lorentz_boost(VX2_IMF, VY2_IMF, VZ2_IMF, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2023 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
2024 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
2025 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
2027 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
2040 lorentz_boost(VX2_IMF, VY2_IMF, VZ2_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
2041 lorentz_boost(vx2ev_imf, vy2ev_imf, vz2ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2042 VIMF_CM[0] = VX2OUT;
2043 VIMF_CM[1] = VY2OUT;
2044 VIMF_CM[2] = VZ2OUT;
2049 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
2050 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2051 VFP1_CM[0] = VX2OUT;
2052 VFP1_CM[1] = VY2OUT;
2053 VFP1_CM[2] = VZ2OUT;
2055 if (FF11 == 0 && FIMF11 == 0)
2068 for (
G4int I = 0; I < 3; I++)
2071 else if (FF11 == 1 && FIMF11 == 0)
2078 opt->optimfallowed = 0;
2087 G4int IEV_TAB_FIS = 0, imode = 0;
2089 G4double vx1_fission = 0., vy1_fission = 0., vz1_fission = 0.;
2090 G4double vx2_fission = 0., vy2_fission = 0., vz2_fission = 0.;
2091 G4double vx_eva_sc = 0., vy_eva_sc = 0., vz_eva_sc = 0.;
2117 for (
int IJ = 0; IJ < IEV_TAB_FIS; IJ++)
2119 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
2120 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
2121 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
2136 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
2137 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
2138 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
2140 IEV_TAB = IEV_TAB + IEV_TAB_FIS;
2152 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
2153 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2154 lorentz_boost(vx1_fission, vy1_fission, vz1_fission, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT, &VZOUT);
2155 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2156 VFP1_CM[0] = VX2OUT;
2157 VFP1_CM[1] = VY2OUT;
2158 VFP1_CM[2] = VZ2OUT;
2167 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
2168 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2169 lorentz_boost(vx2_fission, vy2_fission, vz2_fission, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT, &VZOUT);
2170 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2171 VFP2_CM[0] = VX2OUT;
2172 VFP2_CM[1] = VY2OUT;
2173 VFP2_CM[2] = VZ2OUT;
2175 else if (FF11 == 0 && FIMF11 == 1)
2179 opt->optimfallowed = 0;
2193 G4int NbLamimf1 = 0;
2194 G4double pbH1 = (af - zf) / (af - zf + aimf - zimf);
2195 for (
G4int i = 0; i < NbLamH; i++)
2208 EkinR1 = tkeimf * aimf / (af + aimf);
2209 EkinR2 = tkeimf * af / (af + aimf);
2210 V1 = std::sqrt(EkinR1 / af) * 1.3887;
2211 V2 = std::sqrt(EkinR2 / aimf) * 1.3887;
2213 VPERP1 = std::sqrt(
V1 *
V1 - VZ1_IMFS * VZ1_IMFS);
2215 G4double VX1_IMFS = VPERP1 * std::sin(ALPHA1);
2216 G4double VY1_IMFS = VPERP1 * std::cos(ALPHA1);
2221 EEIMFP = ee * af / (af + aimf);
2222 EEIMF = ee * aimf / (af + aimf);
2225 IINERTTOT = 0.40 * 931.490 * 1.160 * 1.160 * (std::pow(aimf, 5.0 / 3.0) + std::pow(af, 5.0 / 3.0)) +
2226 931.490 * 1.160 * 1.160 * aimf * af / (aimf + af) *
2227 (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.)) *
2228 (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.));
2230 JPRFHEAVY = jprf1 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(af, 5.0 / 3.0) / IINERTTOT;
2231 JPRFLIGHT = jprf1 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(aimf, 5.0 / 3.0) / IINERTTOT;
2233 G4double zffs = 0., affs = 0., vx1ev_imfs = 0., vy1ev_imfs = 0., vz1ev_imfs = 0., jprf3 = 0.;
2257 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
2259 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
2260 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
2261 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
2276 lorentz_boost(vx1ev_imfs, vy1ev_imfs, vz1ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2277 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
2278 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
2279 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
2281 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
2284 opt->optimfallowed = 0;
2290 G4double zffimfs = 0., affimfs = 0., vx2ev_imfs = 0., vy2ev_imfs = 0., vz2ev_imfs = 0., jprf4 = 0.;
2314 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
2316 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
2317 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
2318 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
2333 lorentz_boost(vx2ev_imfs, vy2ev_imfs, vz2ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2334 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
2335 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
2336 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
2338 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
2352 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
2353 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2354 lorentz_boost(VX1_IMFS, VY1_IMFS, VZ1_IMFS, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT, &VZOUT);
2355 lorentz_boost(vx1ev_imfs, vy1ev_imfs, vz1ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2356 VFP1_CM[0] = VX2OUT;
2357 VFP1_CM[1] = VY2OUT;
2358 VFP1_CM[2] = VZ2OUT;
2365 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
2366 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2367 lorentz_boost(VX2_IMFS, VY2_IMFS, VZ2_IMFS, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT, &VZOUT);
2368 lorentz_boost(vx2ev_imfs, vy2ev_imfs, vz2ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2369 VFP2_CM[0] = VX2OUT;
2370 VFP2_CM[1] = VY2OUT;
2371 VFP2_CM[2] = VZ2OUT;
2376 if (ftype != 1 && ftype != 21)
2396 if (IOUNSTABLE == 1)
2403 for (
G4int I = 0; I < ILOOP; I++)
2405 for (
G4int IJ = 0; IJ < 5; IJ++)
2406 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2408 IEV_TAB = IEV_TAB + ILOOP;
2429 if (IOUNSTABLE == 1)
2436 for (
G4int I = 0; I < ILOOP; I++)
2438 for (
G4int IJ = 0; IJ < 5; IJ++)
2439 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2441 IEV_TAB = IEV_TAB + ILOOP;
2462 if (IOUNSTABLE == 1)
2469 for (
G4int I = 0; I < ILOOP; I++)
2471 for (
G4int IJ = 0; IJ < 5; IJ++)
2472 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2474 IEV_TAB = IEV_TAB + ILOOP;
2481 if (ftype == 1 || ftype == 21)
2500 if (IOUNSTABLE == 1)
2507 for (
G4int I = 0; I < ILOOP; I++)
2509 for (
G4int IJ = 0; IJ < 5; IJ++)
2510 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2512 IEV_TAB = IEV_TAB + ILOOP;
2531 if (IOUNSTABLE == 1)
2538 for (
G4int I = 0; I < ILOOP; I++)
2540 for (
G4int IJ = 0; IJ < 5; IJ++)
2541 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2543 IEV_TAB = IEV_TAB + ILOOP;
2564 if (IOUNSTABLE == 1)
2571 for (
G4int I = 0; I < ILOOP; I++)
2573 for (
G4int IJ = 0; IJ < 5; IJ++)
2574 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2576 IEV_TAB = IEV_TAB + ILOOP;
2582 if ((ftype == 1 || ftype == 21) && (AFP2 <= 0 || AFP1 <= 0 || ZFP2 <= 0 || ZFP1 <= 0))
2584 std::cout <<
"ZFP1:" << ZFP1 << std::endl;
2585 std::cout <<
"AFP1:" << AFP1 << std::endl;
2586 std::cout <<
"ZFP2:" << ZFP2 << std::endl;
2587 std::cout <<
"AFP2:" << AFP2 << std::endl;
2591 EV_TAB[IEV_TAB][0] = ZFP1;
2592 EV_TAB[IEV_TAB][1] = AFP1;
2593 EV_TAB[IEV_TAB][5] = SFP1;
2594 EV_TAB[IEV_TAB][2] = VFP1_CM[0];
2595 EV_TAB[IEV_TAB][3] = VFP1_CM[1];
2596 EV_TAB[IEV_TAB][4] = VFP1_CM[2];
2597 IEV_TAB = IEV_TAB + 1;
2601 EV_TAB[IEV_TAB][0] = ZFP2;
2602 EV_TAB[IEV_TAB][1] = AFP2;
2603 EV_TAB[IEV_TAB][5] = SFP2;
2604 EV_TAB[IEV_TAB][2] = VFP2_CM[0];
2605 EV_TAB[IEV_TAB][3] = VFP2_CM[1];
2606 EV_TAB[IEV_TAB][4] = VFP2_CM[2];
2607 IEV_TAB = IEV_TAB + 1;
2612 EV_TAB[IEV_TAB][0] = ZFPIMF;
2613 EV_TAB[IEV_TAB][1] = AFPIMF;
2614 EV_TAB[IEV_TAB][5] = SFPIMF;
2615 EV_TAB[IEV_TAB][2] = VIMF_CM[0];
2616 EV_TAB[IEV_TAB][3] = VIMF_CM[1];
2617 EV_TAB[IEV_TAB][4] = VIMF_CM[2];
2618 IEV_TAB = IEV_TAB + 1;
3886 G4double problamb0 = (*problamb0_par);
4038 G4int choice_fisspart = 0;
4116 G4int fimf_allowed = opt->optimfallowed;
4138 mglms(a, zprf, fiss->optshp, &maz);
4139 mglms(a - 1.0, zprf, fiss->optshp, &ma1z);
4140 mglms(a - 2.0, zprf, fiss->optshp, &ma2z);
4141 mglms(a - 1.0, zprf - 1.0, fiss->optshp, &ma1z1);
4142 mglms(a - 2.0, zprf - 1.0, fiss->optshp, &ma2z1);
4143 mglms(a - 3.0, zprf - 1.0, fiss->optshp, &ma3z1);
4144 mglms(a - 3.0, zprf - 2.0, fiss->optshp, &ma3z2);
4145 mglms(a - 4.0, zprf - 2.0, fiss->optshp, &ma4z2);
4149 mglw(a, zprf, &maz);
4150 mglw(a - 1.0, zprf, &ma1z);
4151 mglw(a - 1.0, zprf - 1.0, &ma1z1);
4152 mglw(a - 2.0, zprf - 1.0, &ma2z1);
4153 mglw(a - 3.0, zprf - 1.0, &ma3z1);
4154 mglw(a - 3.0, zprf - 2.0, &ma3z2);
4155 mglw(a - 4.0, zprf - 2.0, &ma4z2);
4158 if ((a - 1.) == 3.0 && (zprf - 1.0) == 2.0)
4160 if ((a - 1.) == 4.0 && (zprf - 1.0) == 2.0)
4166 sd = ma2z1 - maz - 2.2246;
4167 st = ma3z1 - maz - 8.481977;
4168 she = ma3z2 - maz - 7.7181660;
4169 sa = ma4z2 - maz - 28.295992;
4201 if (zprf <= 1.0e0 || a <= 1.0e0 || (a - zprf) < 0.0)
4214 if (zprf <= 1.0e0 || a <= 2.0e0 || (a - zprf) < 1.0)
4227 if (zprf <= 1.0e0 || a <= 3.0e0 || (a - zprf) < 2.0)
4240 if (a - 4.0 <= 0.0 || zprf <= 2.0 || (a - zprf) < 2.0)
4253 if (a - 3.0 <= 0.0 || zprf <= 2.0 || (a - zprf) < 1.0)
4261 bhe =
max(bhe, 0.1);
4308 barfit(k, k + j, il, &sbfis, &segs, &selmax);
4309 if ((fiss->optshp == 1) || (fiss->optshp == 3))
4311 ef =
G4double(sbfis) - ecld->ecgnz[j][k];
4317 ef = ef * (4.5114 - 2.2687 * (a - zprf) / zprf);
4321 ef = ef * (3.3931 - 1.5338 * (a - zprf) / zprf);
4326 if ((a - zprf) / zprf > 1.52)
4327 ef = ef * (1.1222 - 0.10886 * (a - zprf) / zprf) - 0.1;
4329 if (k >= 94 && k <= 98 && j < 158)
4332 if (
mod(j, 2) == 0 &&
mod(k, 2) == 0)
4336 ef = ef - (11.54108 * (a - zprf) / zprf - 18.074);
4340 if (
mod(j, 2) == 1 &&
mod(k, 2) == 1)
4344 ef = ef - (14.567 * (a - zprf) / zprf - 23.266);
4348 if (
mod(j, 2) == 0 &&
mod(k, 2) == 1)
4352 ef = ef - (13.662 * (a - zprf) / zprf - 21.656);
4356 if (
mod(j, 2) == 1 &&
mod(k, 2) == 0)
4360 ef = ef - (13.662 * (a - zprf) / zprf - 21.656);
4381 ef = ef + 0.51 * (1115. - 938. + sn - slamb0) / std::pow(a, 2. / 3.);
4390 xx =
fissility((k + j), k, NbLam0, sn, slamb0, fiss->optxfis);
4414 bshell = ecld->ecgnz[in][iz] - ecld->vgsld[in][iz];
4415 defbet = ecld->beta2[in][iz];
4417 iinert = 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a, 5.0 / 3.0) * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
4418 erot = jprf * jprf * 197.328 * 197.328 / (2. * iinert);
4421 bsbkbc(a, zprf, &bscn, &bkcn, &bccn);
4425 a, zprf, ee, 0.0, &densg, bshell, bscn, bkcn, &temp, fiss->optshp, fiss->optcol, defbet, &ecor, jprf, 0, &qrcn);
4479 lorb(a, a - 1., jprf, ee - sn, &dlout, &sdlout);
4480 djprf =
gausshaz(1, dlout, sdlout);
4483 jprfn = jprf + djprf;
4484 jprfn =
dint(std::abs(jprfn));
4486 bshell = ecld->ecgnz[ind][izd] - ecld->vgsld[ind][izd];
4487 defbet = ecld->beta2[ind][izd];
4490 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a - 1., 5.0 / 3.0) * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
4491 erotn = jprfn * jprfn * 197.328 * 197.328 / (2. * iinert);
4492 bsbkbc(a - 1., zprf, &bs, &bk, &bc);
4524 std::cout <<
"WARNING: FVMAXHAZ_NEUT CALLED MORE THAN 100 TIMES" << std::endl;
4527 if (ecn > (ee - sn))
4529 if ((ee - sn) < rnt)
4561 lorb(a, a - 1., jprf, ee - sbp, &dlout, &sdlout);
4562 djprf =
gausshaz(1, dlout, sdlout);
4565 jprfp = jprf + djprf;
4566 jprfp =
dint(std::abs(jprfp));
4568 bshell = ecld->ecgnz[ind][izd] - ecld->vgsld[ind][izd];
4569 defbet = ecld->beta2[ind][izd];
4572 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a - 1., 5.0 / 3.0) * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
4573 erotp = jprfp * jprfp * 197.328 * 197.328 / (2. * iinert);
4575 bsbkbc(a - 1., zprf - 1., &bs, &bk, &bc);
4607 std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
4610 if (ecp > (ee - sbp))
4612 if ((ee - sbp) < rpt)
4623 ecp = 2.0 * pt + bp;
4639 if ((in >= 2) && (iz >= 2))
4645 lorb(a, a - 2., jprf, ee - sbd, &dlout, &sdlout);
4646 djprf =
gausshaz(1, dlout, sdlout);
4649 jprfd = jprf + djprf;
4650 jprfd =
dint(std::abs(jprfd));
4652 bshell = ecld->ecgnz[ind][izd] - ecld->vgsld[ind][izd];
4653 defbet = ecld->beta2[ind][izd];
4656 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a - 2., 5.0 / 3.0) * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
4657 erotd = jprfd * jprfd * 197.328 * 197.328 / (2. * iinert);
4659 bsbkbc(a - 2., zprf - 1., &bs, &bk, &bc);
4692 std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
4695 if (ecd > (ee - sbd))
4697 if ((ee - sbd) < rdt)
4708 ecd = 2.0 * dt + bd;
4724 if ((in >= 3) && (iz >= 2))
4730 lorb(a, a - 3., jprf, ee - sbt, &dlout, &sdlout);
4731 djprf =
gausshaz(1, dlout, sdlout);
4734 jprft = jprf + djprf;
4735 jprft =
dint(std::abs(jprft));
4737 bshell = ecld->ecgnz[ind][izd] - ecld->vgsld[ind][izd];
4738 defbet = ecld->beta2[ind][izd];
4741 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a - 3., 5.0 / 3.0) * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
4742 erott = jprft * jprft * 197.328 * 197.328 / (2. * iinert);
4744 bsbkbc(a - 3., zprf - 1., &bs, &bk, &bc);
4777 std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
4780 if (ect > (ee - sbt))
4782 if ((ee - sbt) < rtt)
4793 ect = 2.0 * tt + bt;
4809 if ((in >= 3) && (iz >= 3))
4815 lorb(a, a - 4., jprf, ee - sba, &dlout, &sdlout);
4816 djprf =
gausshaz(1, dlout, sdlout);
4819 jprfa = jprf + djprf;
4820 jprfa =
dint(std::abs(jprfa));
4822 bshell = ecld->ecgnz[ind][izd] - ecld->vgsld[ind][izd];
4823 defbet = ecld->beta2[ind][izd];
4826 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a - 4., 5.0 / 3.0) * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
4827 erota = jprfa * jprfa * 197.328 * 197.328 / (2. * iinert);
4829 bsbkbc(a - 4., zprf - 2., &bs, &bk, &bc);
4862 std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
4865 if (eca > (ee - sba))
4867 if ((ee - sba) < rat)
4878 eca = 2.0 * at + ba;
4894 if ((in >= 2) && (iz >= 3))
4900 lorb(a, a - 3., jprf, ee - sbhe, &dlout, &sdlout);
4901 djprf =
gausshaz(1, dlout, sdlout);
4904 jprfhe = jprf + djprf;
4905 jprfhe =
dint(std::abs(jprfhe));
4907 bshell = ecld->ecgnz[ind][izd] - ecld->vgsld[ind][izd];
4908 defbet = ecld->beta2[ind][izd];
4911 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a - 3., 5.0 / 3.0) * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
4912 erothe = jprfhe * jprfhe * 197.328 * 197.328 / (2. * iinert);
4914 bsbkbc(a - 3., zprf - 2., &bs, &bk, &bc);
4947 std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
4950 if (eche > (ee - sbhe))
4952 if ((ee - sbhe) < rhet)
4963 eche = 2.0 * het + bhe;
4980 if (in >= 2 && NbLam0 > 0)
4986 lorb(a, a - 1., jprf, ee - slamb0, &dlout, &sdlout);
4987 djprf =
gausshaz(1, dlout, sdlout);
4990 jprflamb0 = jprf + djprf;
4991 jprflamb0 =
dint(std::abs(jprflamb0));
4993 bshell = ecld->ecgnz[ind][izd] - ecld->vgsld[ind][izd];
4994 defbet = ecld->beta2[ind][izd];
4997 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a - 1., 5.0 / 3.0) * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
4998 erotlamb0 = jprflamb0 * jprflamb0 * 197.328 * 197.328 / (2. * iinert);
4999 bsbkbc(a - 1., zprf, &bs, &bk, &bc);
5031 std::cout <<
"WARNING: FVMAXHAZ_NEUT CALLED MORE THAN 100 TIMES" << std::endl;
5034 if (eclamb0 > (ee - slamb0))
5036 if ((ee - slamb0) < rlamb0t)
5037 eclamb0 = ee - slamb0;
5046 eclamb0 = 2.0 * lamb0t;
5078 gn =
width(a, zprf, 1.0, 0.0, nt, 0.0, sn, ee - erotn) * densn / densg;
5086 gp =
width(a, zprf, 1.0, 1.0, pt, bp, sbp, ee - erotp) * densp / densg *
pen(a, 1.0, omegap, pt);
5094 gd =
width(a, zprf, 2.0, 1.0, dt, bd, sbd, ee - erotd) * densd / densg *
pen(a, 2.0, omegad, dt);
5102 gt =
width(a, zprf, 3.0, 1.0, tt, bt, sbt, ee - erott) * denst / densg *
pen(a, 3.0, omegat, tt);
5110 ghe =
width(a, zprf, 3.0, 2.0, het, bhe, sbhe, ee - erothe) * denshe / densg *
pen(a, 3.0, omegahe, het);
5118 ga =
width(a, zprf, 4.0, 2.0, at, ba, sba, ee - erota) * densa / densg *
pen(a, 4.0, omegaa, at);
5120 if (denslamb0 <= 0.0)
5126 glamb0 =
width(a, zprf, 1.0, -2.0, lamb0t, 0.0, slamb0, ee - erotlamb0) * denslamb0 / densg;
5134 G4int izcn = 0, incn = 0, inmin = 0, inmax = 0, inmi = 0, inma = 0;
5137 if (fimf_allowed == 0 || zprf <= 5.0 || a <= 7.0)
5146 mglms(a, zprf, opt->optshpimf, &mazz);
5161 inmin =
max(inmin, incn - inma);
5162 inmax =
min(inmax, incn - inmi);
5164 inmax =
max(inmax, inmin);
5166 for (
G4int iaimf = izimf + inmin; iaimf <= izimf + inmax; iaimf++)
5169 if (aimf >= a || zimf >= zprf)
5176 mglms(a - aimf, zprf - zimf, opt->optshpimf, &mares);
5177 mglms(aimf, zimf, opt->optshpimf, &maimf);
5182 defbetimf = ecld->beta2[
idnint(aimf - zimf)][
idnint(zimf)] +
5183 ecld->beta2[
idnint(a - aimf - zprf + zimf)][
idnint(zprf - zimf)];
5185 iinert = 0.40 * 931.490 * 1.160 * 1.160 * std::pow(a, 5.0 / 3.0) *
5186 (std::pow(aimf, 5.0 / 3.0) + std::pow(a - aimf, 5.0 / 3.0)) +
5187 931.490 * 1.160 * 1.160 * aimf * (a - aimf) / a *
5188 (std::pow(aimf, 1.0 / 3.0) + std::pow(a - aimf, 1.0 / 3.0)) *
5189 (std::pow(aimf, 1.0 / 3.0) + std::pow(a - aimf, 1.0 / 3.0));
5191 erot = jprf * jprf * 197.328 * 197.328 / (2.0 * iinert);
5194 if (densg == 0.0 || ee < (sbimf + erot))
5205 a, zprf, ee, sbimf, &densimf, 0.0, bsimf, 1.0, &timf, 0, 0, defbetimf, &ecor, jprf, 2, &qr);
5207 imfarg = (sbimf + erotcn - erot) / timf;
5221 width(a, zprf, aimf, zimf, timf, bimf, sbimf, ee - erot) * std::exp(-imfarg) * qr / qrcn;
5224 gimf3 = gimf3 + width_imf;
5241 inmin =
max(inmin, incn - inma);
5242 inmax =
min(inmax, incn - inmi);
5244 inmax =
max(inmax, inmin);
5246 for (
G4int iaimf = izimf + inmin; iaimf <= izimf + inmax; iaimf++)
5249 if (aimf >= a || zimf >= zprf)
5256 mglms(a - aimf, zprf - zimf, opt->optshpimf, &mares);
5257 mglms(aimf, zimf, opt->optshpimf, &maimf);
5262 defbetimf = ecld->beta2[
idnint(aimf - zimf)][
idnint(zimf)] +
5263 ecld->beta2[
idnint(a - aimf - zprf + zimf)][
idnint(zprf - zimf)];
5265 iinert = 0.40 * 931.490 * 1.160 * 1.160 * std::pow(a, 5.0 / 3.0) *
5266 (std::pow(aimf, 5.0 / 3.0) + std::pow(a - aimf, 5.0 / 3.0)) +
5267 931.490 * 1.160 * 1.160 * aimf * (a - aimf) / a *
5268 (std::pow(aimf, 1.0 / 3.0) + std::pow(a - aimf, 1.0 / 3.0)) *
5269 (std::pow(aimf, 1.0 / 3.0) + std::pow(a - aimf, 1.0 / 3.0));
5271 erot = jprf * jprf * 197.328 * 197.328 / (2.0 * iinert);
5274 if (densg == 0.0 || ee < (sbimf + erot))
5285 a, zprf, ee, sbimf, &densimf, 0.0, bsimf, 1.0, &timf, 0, 0, defbetimf, &ecor, jprf, 2, &qr);
5287 imfarg = (sbimf + erotcn - erot) / timf;
5299 width_imf =
width(a, zprf, aimf, zimf, timf, bimf, sbimf, ee - erot) * std::exp(-imfarg) * qr /
5303 gimf5 = gimf5 + width_imf;
5308 if (gimf3 <= 0.0 || gimf5 <= 0.0)
5317 b_imf = (std::log10(gimf3) - std::log10(gimf5)) / (std::log10(3.0) - std::log10(5.0));
5321 if (b_imf <= -100.0)
5329 a_imf = gimf3 / std::pow(3.0, b_imf);
5330 gimf = a_imf * (std::pow(zprf, b_imf + 1.0) - std::pow(3.0, b_imf + 1.0)) / (b_imf + 1.0);
5342 pa = (ald->av) * a + (ald->as) * std::pow(a, 2. / 3.) + (ald->ak) * std::pow(a, 1. / 3.);
5343 gamma = 2.5 * pa * std::pow(a, -4. / 3.);
5344 gfactor = 1. + gamma * ecld->ecgnz[in][iz];
5350 gtemp = 17.60 / (std::pow(a, 0.699) * std::sqrt(gfactor));
5354 gg = 0.624e-9 * std::pow(a, 1.6) * std::pow(gtemp, 5.);
5358 if (gammaemission == 1)
5365 gsum = ga + ghe + gd + gt + gp + gn + gimf + gg + glamb0;
5380 if (fiss->ifis == 0 || (zprf * zprf / a <= 22.74 && zprf < 60.))
5389 fission_width(zprf, a, ee, bssp, bksp, ef, y, &gf, &temp, jprf, 0, 1, fiss->optcol, fiss->optshp, densg);
5409 if (fiss->bet <= 0.)
5411 gtotal = ga + ghe + gp + gd + gt + gn + gg + gimf + gf + glamb0;
5428 probf = gf / gtotal;
5429 probn = gn / gtotal;
5430 probp = gp / gtotal;
5431 probd = gd / gtotal;
5432 probt = gt / gtotal;
5433 probhe = ghe / gtotal;
5434 proba = ga / gtotal;
5435 probg = gg / gtotal;
5436 probimf = gimf / gtotal;
5437 problamb0 = glamb0 / gtotal;
5453 fomega_sp(a, y, &mfcd, &omegasp, &homegasp);
5454 cf =
cram((NbLam0 > 0 ? fiss->bethyp : fiss->bet), homegasp);
5457 fomega_gs(a, zprf, &k1, &omegags, &homegags);
5458 tauc =
tau((NbLam0 > 0 ? fiss->bethyp : fiss->bet), homegags, ef, ft);
5471 part_fiss((NbLam0 > 0 ? fiss->bethyp : fiss->bet), gsum, gf, y, tauc, ts1, tsum, &choice_fisspart, zprf, a, ft, &t_lapse, &gff);
5476 tsum = tsum + t_lapse;
5479 if (choice_fisspart == 2)
5499 gtotal = ga + ghe + gp + gd + gt + gn + gimf + gg + glamb0;
5517 probn = gn / gtotal;
5518 probp = gp / gtotal;
5519 probd = gd / gtotal;
5520 probt = gt / gtotal;
5521 probhe = ghe / gtotal;
5522 proba = ga / gtotal;
5523 probg = gg / gtotal;
5524 probimf = gimf / gtotal;
5525 problamb0 = glamb0 / gtotal;
5531 gtotal = ga + ghe + gp + gd + gt + gn + gg + gimf + glamb0;
5548 probn = gn / gtotal;
5549 probp = gp / gtotal;
5550 probd = gd / gtotal;
5551 probt = gt / gtotal;
5552 probhe = ghe / gtotal;
5553 proba = ga / gtotal;
5554 probg = gg / gtotal;
5555 probimf = gimf / gtotal;
5556 problamb0 = glamb0 / gtotal;
5560 ptotl = probp + probd + probt + probn + probhe + proba + probg + probimf + probf + problamb0;
5566 (*probp_par) = probp;
5567 (*probd_par) = probd;
5568 (*probt_par) = probt;
5569 (*probn_par) = probn;
5570 (*probhe_par) = probhe;
5571 (*proba_par) = proba;
5572 (*probg_par) = probg;
5573 (*probimf_par) = probimf;
5574 (*problamb0_par) = problamb0;
5575 (*probf_par) = probf;
5576 (*ptotl_par) = ptotl;
5583 (*slamb0_par) = slamb0;
5596 (*eclamb0_par) = eclamb0;
5604 (*jprfn_par) = jprfn;
5605 (*jprfp_par) = jprfp;
5606 (*jprfd_par) = jprfd;
5607 (*jprft_par) = jprft;
5608 (*jprfhe_par) = jprfhe;
5609 (*jprfa_par) = jprfa;
5610 (*jprflamb0_par) = jprflamb0;
6777 for (
G4int init_i = 0; init_i < 7; init_i++)
6782 for (
G4int init_i = 0; init_i < 10; init_i++)
6787 G4double a = 0.0, z = 0.0, amin = 0.0, amax = 0.0, amin2 = 0.0;
6788 G4double amax2 = 0.0, aa = 0.0, zz = 0.0, bfis = 0.0;
6789 G4double bfis0 = 0.0, ell = 0.0, el = 0.0, egs = 0.0, el80 = 0.0, el20 = 0.0;
6790 G4double elmax = 0.0, sel80 = 0.0, sel20 = 0.0, x = 0.0, y = 0.0, q = 0.0, qa = 0.0, qb = 0.0;
6791 G4double aj = 0.0, ak = 0.0, a1 = 0.0, a2 = 0.0;
6793 G4int i = 0, j = 0, k = 0, m = 0;
6796 G4double emncof[4][5] = { { -9.01100e+2, -1.40818e+3, 2.77000e+3, -7.06695e+2, 8.89867e+2 },
6797 { 1.35355e+4, -2.03847e+4, 1.09384e+4, -4.86297e+3, -6.18603e+2 },
6798 { -3.26367e+3, 1.62447e+3, 1.36856e+3, 1.31731e+3, 1.53372e+2 },
6799 { 7.48863e+3, -1.21581e+4, 5.50281e+3, -1.33630e+3, 5.05367e-2 } };
6801 G4double elmcof[4][5] = { { 1.84542e+3, -5.64002e+3, 5.66730e+3, -3.15150e+3, 9.54160e+2 },
6802 { -2.24577e+3, 8.56133e+3, -9.67348e+3, 5.81744e+3, -1.86997e+3 },
6803 { 2.79772e+3, -8.73073e+3, 9.19706e+3, -4.91900e+3, 1.37283e+3 },
6804 { -3.01866e+1, 1.41161e+3, -2.85919e+3, 2.13016e+3, -6.49072e+2 } };
6806 G4double emxcof[4][6] = { { 9.43596e4, -2.241997e5, 2.223237e5, -1.324408e5, 4.68922e4, -8.83568e3 },
6807 { -1.655827e5, 4.062365e5, -4.236128e5, 2.66837e5, -9.93242e4, 1.90644e4 },
6808 { 1.705447e5, -4.032e5, 3.970312e5, -2.313704e5, 7.81147e4, -1.322775e4 },
6809 { -9.274555e4, 2.278093e5, -2.422225e5, 1.55431e5, -5.78742e4, 9.97505e3 } };
6812 { 5.11819909e+5, -1.30303186e+6, 1.90119870e+6, -1.20628242e+6, 5.68208488e+5, 5.48346483e+4, -2.45883052e+4 },
6813 { -1.13269453e+6, 2.97764590e+6, -4.54326326e+6, 3.00464870e+6, -1.44989274e+6, -1.02026610e+5, 6.27959815e+4 },
6814 { 1.37543304e+6, -3.65808988e+6, 5.47798999e+6, -3.78109283e+6, 1.84131765e+6, 1.53669695e+4, -6.96817834e+4 },
6815 { -8.56559835e+5, 2.48872266e+6, -4.07349128e+6, 3.12835899e+6, -1.62394090e+6, 1.19797378e+5, 4.25737058e+4 },
6816 { 3.28723311e+5, -1.09892175e+6, 2.03997269e+6, -1.77185718e+6, 9.96051545e+5, -1.53305699e+5, -1.12982954e+4 },
6817 { 4.15850238e+4, 7.29653408e+4, -4.93776346e+5, 6.01254680e+5, -4.01308292e+5, 9.65968391e+4, -3.49596027e+3 },
6818 { -1.82751044e+5, 3.91386300e+5, -3.03639248e+5, 1.15782417e+5, -4.24399280e+3, -6.11477247e+3, 3.66982647e+2 }
6821 const G4int sizex = 5;
6822 const G4int sizey = 6;
6823 const G4int sizez = 4;
6825 G4double egscof[sizey][sizey][sizez];
6827 G4double egs1[sizey][sizex] = { { 1.927813e5, 7.666859e5, 6.628436e5, 1.586504e5, -7.786476e3 },
6828 { -4.499687e5, -1.784644e6, -1.546968e6, -4.020658e5, -3.929522e3 },
6829 { 4.667741e5, 1.849838e6, 1.641313e6, 5.229787e5, 5.928137e4 },
6830 { -3.017927e5, -1.206483e6, -1.124685e6, -4.478641e5, -8.682323e4 },
6831 { 1.226517e5, 5.015667e5, 5.032605e5, 2.404477e5, 5.603301e4 },
6832 { -1.752824e4, -7.411621e4, -7.989019e4, -4.175486e4, -1.024194e4 } };
6834 G4double egs2[sizey][sizex] = { { -6.459162e5, -2.903581e6, -3.048551e6, -1.004411e6, -6.558220e4 },
6835 { 1.469853e6, 6.564615e6, 6.843078e6, 2.280839e6, 1.802023e5 },
6836 { -1.435116e6, -6.322470e6, -6.531834e6, -2.298744e6, -2.639612e5 },
6837 { 8.665296e5, 3.769159e6, 3.899685e6, 1.520520e6, 2.498728e5 },
6838 { -3.302885e5, -1.429313e6, -1.512075e6, -6.744828e5, -1.398771e5 },
6839 { 4.958167e4, 2.178202e5, 2.400617e5, 1.167815e5, 2.663901e4 } };
6841 G4double egs3[sizey][sizex] = { { 3.117030e5, 1.195474e6, 9.036289e5, 6.876190e4, -6.814556e4 },
6842 { -7.394913e5, -2.826468e6, -2.152757e6, -2.459553e5, 1.101414e5 },
6843 { 7.918994e5, 3.030439e6, 2.412611e6, 5.228065e5, 8.542465e3 },
6844 { -5.421004e5, -2.102672e6, -1.813959e6, -6.251700e5, -1.184348e5 },
6845 { 2.370771e5, 9.459043e5, 9.026235e5, 4.116799e5, 1.001348e5 },
6846 { -4.227664e4, -1.738756e5, -1.795906e5, -9.292141e4, -2.397528e4 } };
6848 G4double egs4[sizey][sizex] = { { -1.072763e5, -5.973532e5, -6.151814e5, 7.371898e4, 1.255490e5 },
6849 { 2.298769e5, 1.265001e6, 1.252798e6, -2.306276e5, -2.845824e5 },
6850 { -2.093664e5, -1.100874e6, -1.009313e6, 2.705945e5, 2.506562e5 },
6851 { 1.274613e5, 6.190307e5, 5.262822e5, -1.336039e5, -1.115865e5 },
6852 { -5.715764e4, -2.560989e5, -2.228781e5, -3.222789e3, 1.575670e4 },
6853 { 1.189447e4, 5.161815e4, 4.870290e4, 1.266808e4, 2.069603e3 } };
6855 for (i = 0; i < sizey; i++)
6857 for (j = 0; j < sizex; j++)
6859 egscof[i][j][0] = egs1[i][j];
6860 egscof[i][j][1] = egs2[i][j];
6861 egscof[i][j][2] = egs3[i][j];
6862 egscof[i][j][3] = egs4[i][j];
6867 if (iz < 19 || iz > 122)
6872 if (iz > 122 && il > 0)
6880 amin = 1.2e0 * z + 0.01e0 * z * z;
6881 amax = 5.8e0 * z - 0.024e0 * z * z;
6883 if (a < amin || a > amax)
6896 for (i = 0; i < 7; i++)
6898 for (j = 0; j < 7; j++)
6900 bfis0 = bfis0 + elzcof[j][i] * pz[i] * pa[j];
6912 amin2 = 1.4e0 * z + 0.009e0 * z * z;
6913 amax2 = 20.e0 + 3.0e0 * z;
6915 if ((a < amin2 - 5.e0 || a > amax2 + 10.e0) && il > 0)
6926 for (i = 0; i < 4; i++)
6928 for (j = 0; j < 5; j++)
6930 el80 = el80 + elmcof[i][j] * pz[j] * pa[i];
6931 el20 = el20 + emncof[i][j] * pz[j] * pa[i];
6942 for (i = 0; i < 4; i++)
6944 for (j = 0; j < 6; j++)
6946 elmax = elmax + emxcof[i][j] * pz[j] * pa[i];
6958 x = sel20 / (*selmax);
6959 y = sel80 / (*selmax);
6964 q = 0.2 / (std::pow(sel20, 2) * std::pow(sel80, 2) * (sel20 - sel80));
6965 qa = q * (4.0 * std::pow(sel80, 3) - std::pow(sel20, 3));
6966 qb = -q * (4.0 * std::pow(sel80, 2) - std::pow(sel20, 2));
6967 bfis = bfis * (1.0 + qa * std::pow(el, 2) + qb * std::pow(el, 3));
6972 aj = (-20.0 * std::pow(x, 5) + 25.e0 * std::pow(x, 4) - 4.0) * std::pow((y - 1.0), 2) * y * y;
6973 ak = (-20.0 * std::pow(y, 5) + 25.0 * std::pow(y, 4) - 1.0) * std::pow((x - 1.0), 2) * x * x;
6974 q = 0.2 / (std::pow((y - x) * ((1.0 - x) * (1.0 - y) * x * y), 2));
6975 qa = q * (aj * y - ak * x);
6976 qb = -q * (aj * (2.0 * y + 1.0) - ak * (2.0 * x + 1.0));
6978 a1 = 4.0 * std::pow(z, 5) - 5.0 * std::pow(z, 4) + 1.0;
6979 a2 = qa * (2.e0 * z + 1.e0);
6980 bfis = bfis * (a1 + (z - 1.e0) * (a2 + qb * z) * z * z * (z - 1.e0));
7000 for (k = 0; k < 4; k++)
7002 for (l = 0; l < 6; l++)
7004 for (m = 0; m < 5; m++)
7006 egs = egs + egscof[l][m][k] * pz[l] * pa[k] * pl[2 * m];
9799 Delta_U1_shell_max = -2.45;
9805 Delta_U2_shell = -2.45;
9815 G4double Fwidth_asymm1, Fwidth_asymm2, Fwidth_symm;
9817 Fwidth_asymm1 = 0.65;
9818 Fwidth_asymm2 = 0.65;
9860 Friction_factor = 1.0;
9863 G4double cN_asymm1_shell, cN_asymm2_shell;
9864 G4double gamma, gamma_heavy1, gamma_heavy2;
9880 G4double Sasymm1 = 0., Sasymm2 = 0., Ssymm = 0., Ysum = 0., Yasymm = 0.;
9882 G4double wNasymm1_saddle, wNasymm2_saddle, wNsymm_saddle;
9883 G4double wNasymm2_scission, wNsymm_scission;
9884 G4double wNasymm1, wNasymm2, wNsymm;
9898 G4double beta = 0., beta1 = 0., beta2 = 0.;
9907 G4double A_levdens_heavy1, A_levdens_heavy2;
9911 G4double epsilon_1_saddle, epsilon0_1_saddle;
9912 G4double epsilon_2_saddle, epsilon0_2_saddle, epsilon_symm_saddle;
9917 G4double E_eff1_saddle, E_eff2_saddle;
9918 G4double Epot0_mode1_saddle, Epot0_mode2_saddle, Epot0_symm_saddle;
9919 G4double Epot_mode1_saddle, Epot_mode2_saddle, Epot_symm_saddle;
9920 G4double E_defo, E_defo1, E_defo2, E_scission_pre = 0., E_scission_post;
9926 G4double MassCurv_scis, MassCurv_sadd;
9928 G4double Nheavy1_shell, Nheavy2_shell;
9933 G4double Z_scission, N_scission, A_scission;
9935 G4double beta1gs = 0., beta2gs = 0., betags = 0.;
9937 G4double DSN132, Delta_U1_shell, E_eff0_saddle;
9938 G4int NbLam0 = (*NbLam0_par);
9943 cN_asymm1_shell = 0.700 *
N / Z;
9944 cN_asymm2_shell = 0.040 *
N / Z;
9948 DSN132 = Nheavy1_in -
N / Z * Zheavy1_in;
9949 Aheavy1 = Nheavy1_in + Zheavy1_in + 0.340 * DSN132;
9958 Delta_U1_shell = Delta_U1_shell_max + U1NZ_SLOPE * std::abs(DSN132);
9959 Delta_U1_shell =
min(0., Delta_U1_shell);
9963 Nheavy1 =
N /
A * Aheavy1;
9964 Aheavy2 = Nheavy2 *
A /
N;
9968 A_levdens =
A / xLevdens;
9969 gamma = A_levdens / (0.40 * std::pow(
A, 1.3333)) * FGAMMA;
9970 A_levdens_heavy1 = Aheavy1 / xLevdens;
9971 gamma_heavy1 = A_levdens_heavy1 / (0.40 * std::pow(Aheavy1, 1.3333)) * FGAMMA * FGAMMA1;
9972 A_levdens_heavy2 = Aheavy2 / xLevdens;
9973 gamma_heavy2 = A_levdens_heavy2 / (0.40 * std::pow(Aheavy2, 1.3333)) * FGAMMA;
9977 E_saddle_scission = (-24. + 0.02227 * Z * Z / std::pow(
A, 0.33333)) * Friction_factor;
9978 E_saddle_scission =
max(0.0, E_saddle_scission);
9984 Z2_over_A_eff = Z * Z /
A;
9986 if (Z2_over_A_eff < 34.0)
9987 MassCurv_scis = std::pow(10., -1.093364 + 0.082933 * Z2_over_A_eff - 0.0002602 * Z2_over_A_eff * Z2_over_A_eff);
9989 MassCurv_scis = std::pow(10., 3.053536 - 0.056477 * Z2_over_A_eff + 0.0002454 * Z2_over_A_eff * Z2_over_A_eff);
9994 MassCurv_sadd = X_s2s * MassCurv_scis;
9996 cN_symm = 8.0 / std::pow(
N, 2.) * MassCurv_scis;
9997 cN_symm_sadd = 8.0 / std::pow(
N, 2.) * MassCurv_sadd;
9999 Nheavy1_shell = Nheavy1;
10002 Nheavy1_eff = (cN_symm_sadd * Nsymm +
10003 cN_asymm1_shell *
Uwash(E /
A * Aheavy1, Ecrit, FREDSHELL, gamma_heavy1) * Nheavy1_shell) /
10004 (cN_symm_sadd + cN_asymm1_shell *
Uwash(E /
A * Aheavy1, Ecrit, FREDSHELL, gamma_heavy1));
10006 Nheavy1_eff = (cN_symm_sadd * Nsymm + cN_asymm1_shell * Nheavy1_shell) / (cN_symm_sadd + cN_asymm1_shell);
10009 Nheavy2_NZ = Nheavy2;
10010 Nheavy2_shell = Nheavy2_NZ;
10012 Nheavy2_eff = (cN_symm_sadd * Nsymm +
10013 cN_asymm2_shell *
Uwash(E /
A * Aheavy2, Ecrit, FREDSHELL, gamma_heavy2) * Nheavy2_shell) /
10014 (cN_symm_sadd + cN_asymm2_shell *
Uwash(E /
A * Aheavy2, Ecrit, FREDSHELL, gamma_heavy2));
10016 Nheavy2_eff = (cN_symm_sadd * Nsymm + cN_asymm2_shell * Nheavy2_shell) / (cN_symm_sadd + cN_asymm2_shell);
10018 Delta_U1 = Delta_U1_shell + (Nheavy1_shell - Nheavy1_eff) * (Nheavy1_shell - Nheavy1_eff) *
10020 Delta_U1 =
min(Delta_U1, 0.0);
10021 Delta_U2 = Delta_U2_shell + (Nheavy2_shell - Nheavy2_eff) * (Nheavy2_shell - Nheavy2_eff) *
10023 Delta_U2 =
min(Delta_U2, 0.0);
10027 Epot0_mode1_saddle = (Nheavy1_eff - Nsymm) * (Nheavy1_eff - Nsymm) * cN_symm_sadd;
10028 Epot0_mode2_saddle = (Nheavy2_eff - Nsymm) * (Nheavy2_eff - Nsymm) * cN_symm_sadd;
10029 Epot0_symm_saddle = 0.0;
10033 Epot_mode1_saddle = Epot0_mode1_saddle + Delta_U1;
10034 Epot_mode2_saddle = Epot0_mode2_saddle + Delta_U2;
10035 Epot_symm_saddle = Epot0_symm_saddle;
10038 dUeff =
min(Epot_mode1_saddle, Epot_mode2_saddle);
10039 dUeff =
min(dUeff, Epot_symm_saddle);
10040 dUeff = dUeff - Epot_symm_saddle;
10049 epsilon0_1_saddle = Eld - Epot0_mode1_saddle;
10050 epsilon0_2_saddle = Eld - Epot0_mode2_saddle;
10053 epsilon_1_saddle = Eld - Epot_mode1_saddle;
10054 epsilon_2_saddle = Eld - Epot_mode2_saddle;
10056 epsilon_symm_saddle = Eld - Epot_symm_saddle;
10059 eexc1_saddle = epsilon_1_saddle;
10060 eexc2_saddle = epsilon_2_saddle;
10063 EEXC_MAX =
max(eexc1_saddle, eexc2_saddle);
10064 EEXC_MAX =
max(EEXC_MAX, Eld);
10067 epsilon_1_scission = Eld + E_saddle_scission - Epot_mode1_saddle;
10068 epsilon_2_scission = Eld + E_saddle_scission - Epot_mode2_saddle;
10071 epsilon_symm_scission = Eld + E_saddle_scission - Epot_symm_saddle;
10075 epsilon0_1_saddle - Delta_U1 *
Uwash(epsilon_1_saddle /
A * Aheavy1, Ecrit, FREDSHELL, gamma_heavy1);
10077 if (E_eff1_saddle < A_levdens * hbom1 * hbom1)
10078 E_eff1_saddle = A_levdens * hbom1 * hbom1;
10080 wNasymm1_saddle = std::sqrt(
10081 0.50 * std::sqrt(1.0 / A_levdens * E_eff1_saddle) /
10082 (cN_asymm1_shell *
Uwash(epsilon_1_saddle /
A * Aheavy1, Ecrit, FREDSHELL, gamma_heavy1) + cN_symm_sadd));
10085 epsilon0_2_saddle - Delta_U2 *
Uwash(epsilon_2_saddle /
A * Aheavy2, Ecrit, FREDSHELL, gamma_heavy2);
10087 if (E_eff2_saddle < A_levdens * hbom2 * hbom2)
10088 E_eff2_saddle = A_levdens * hbom2 * hbom2;
10090 wNasymm2_saddle = std::sqrt(
10091 0.50 * std::sqrt(1.0 / A_levdens * E_eff2_saddle) /
10092 (cN_asymm2_shell *
Uwash(epsilon_2_saddle /
A * Aheavy2, Ecrit, FREDSHELL, gamma_heavy2) + cN_symm_sadd));
10094 E_eff0_saddle = epsilon_symm_saddle;
10095 if (E_eff0_saddle < A_levdens * hbom3 * hbom3)
10096 E_eff0_saddle = A_levdens * hbom3 * hbom3;
10098 wNsymm_saddle = std::sqrt(0.50 * std::sqrt(1.0 / A_levdens * E_eff0_saddle) / cN_symm_sadd);
10100 if (epsilon_symm_scission > 0.0)
10102 E_HELP =
max(E_saddle_scission, epsilon_symm_scission);
10103 wNsymm_scission = std::sqrt(0.50 * std::sqrt(1.0 / A_levdens * (E_HELP)) / cN_symm);
10107 wNsymm_scission = std::sqrt(0.50 * std::sqrt(1.0 / A_levdens * E_saddle_scission) / cN_symm);
10113 if (E_saddle_scission == 0.0)
10115 wNasymm1_scission = wNasymm1_saddle;
10116 wNasymm2_scission = wNasymm2_saddle;
10120 if (Nheavy1_eff > 75.0)
10122 wNasymm1_scission = std::sqrt(21.0) *
N /
A;
10123 wNasymm2_scission =
max(12.8 - 1.0 * (92.0 - Nheavy2_eff), 1.0) *
N /
A;
10127 wNasymm1_scission = wNasymm1_saddle;
10128 wNasymm2_scission = wNasymm2_saddle;
10132 wNasymm1_scission =
max(wNasymm1_scission, wNasymm1_saddle);
10133 wNasymm2_scission =
max(wNasymm2_scission, wNasymm2_saddle);
10135 wNasymm1 = wNasymm1_scission * Fwidth_asymm1;
10136 wNasymm2 = wNasymm2_scission * Fwidth_asymm2;
10137 wNsymm = wNsymm_scission * Fwidth_symm;
10140 Aheavy1_eff = Nheavy1_eff *
A /
N;
10141 Aheavy2_eff = Nheavy2_eff *
A /
N;
10143 A_levdens_heavy1 = Aheavy1_eff / xLevdens;
10144 A_levdens_heavy2 = Aheavy2_eff / xLevdens;
10145 gamma_heavy1 = A_levdens_heavy1 / (0.40 * std::pow(Aheavy1_eff, 1.3333)) * FGAMMA * FGAMMA1;
10146 gamma_heavy2 = A_levdens_heavy2 / (0.40 * std::pow(Aheavy2_eff, 1.3333)) * FGAMMA;
10148 if (epsilon_symm_saddle < A_levdens * hbom3 * hbom3)
10149 Ssymm = 2.0 * std::sqrt(A_levdens * A_levdens * hbom3 * hbom3) +
10150 (epsilon_symm_saddle - A_levdens * hbom3 * hbom3) / hbom3;
10152 Ssymm = 2.0 * std::sqrt(A_levdens * epsilon_symm_saddle);
10156 if (epsilon0_1_saddle < A_levdens * hbom1 * hbom1)
10157 Ssymm_mode1 = 2.0 * std::sqrt(A_levdens * A_levdens * hbom1 * hbom1) +
10158 (epsilon0_1_saddle - A_levdens * hbom1 * hbom1) / hbom1;
10160 Ssymm_mode1 = 2.0 * std::sqrt(A_levdens * epsilon0_1_saddle);
10162 if (epsilon0_2_saddle < A_levdens * hbom2 * hbom2)
10163 Ssymm_mode2 = 2.0 * std::sqrt(A_levdens * A_levdens * hbom2 * hbom2) +
10164 (epsilon0_2_saddle - A_levdens * hbom2 * hbom2) / hbom2;
10166 Ssymm_mode2 = 2.0 * std::sqrt(A_levdens * epsilon0_2_saddle);
10168 if (epsilon0_1_saddle - Delta_U1 *
Uwash(epsilon_1_saddle /
A * Aheavy1, Ecrit, FREDSHELL, gamma_heavy1) <
10169 A_levdens * hbom1 * hbom1)
10171 2.0 * std::sqrt(A_levdens * A_levdens * hbom1 * hbom1) +
10172 (epsilon0_1_saddle - Delta_U1 *
Uwash(epsilon_1_saddle /
A * Aheavy1, Ecrit, FREDSHELL, gamma_heavy1) -
10173 A_levdens * hbom1 * hbom1) /
10176 Sasymm1 = 2.0 * std::sqrt(A_levdens *
10177 (epsilon0_1_saddle -
10178 Delta_U1 *
Uwash(epsilon_1_saddle /
A * Aheavy1, Ecrit, FREDSHELL, gamma_heavy1)));
10180 if (epsilon0_2_saddle - Delta_U2 *
Uwash(epsilon_2_saddle /
A * Aheavy2, Ecrit, FREDSHELL, gamma_heavy2) <
10181 A_levdens * hbom2 * hbom2)
10183 2.0 * std::sqrt(A_levdens * A_levdens * hbom2 * hbom2) +
10184 (epsilon0_1_saddle - Delta_U1 *
Uwash(epsilon_2_saddle /
A * Aheavy2, Ecrit, FREDSHELL, gamma_heavy2) -
10185 A_levdens * hbom2 * hbom2) /
10188 Sasymm2 = 2.0 * std::sqrt(A_levdens *
10189 (epsilon0_2_saddle -
10190 Delta_U2 *
Uwash(epsilon_2_saddle /
A * Aheavy2, Ecrit, FREDSHELL, gamma_heavy2)));
10192 Yasymm1 = (std::exp(Sasymm1 - Ssymm) - std::exp(Ssymm_mode1 - Ssymm)) * wNasymm1_saddle / wNsymm_saddle * 2.0;
10194 Yasymm2 = (std::exp(Sasymm2 - Ssymm) - std::exp(Ssymm_mode2 - Ssymm)) * wNasymm2_saddle / wNsymm_saddle * 2.0;
10196 Ysum = Ysymm + Yasymm1 + Yasymm2;
10200 Ysymm = Ysymm / Ysum;
10201 Yasymm1 = Yasymm1 / Ysum;
10202 Yasymm2 = Yasymm2 / Ysum;
10203 Yasymm = Yasymm1 + Yasymm2;
10211 if ((epsilon_symm_saddle < epsilon_1_saddle) && (epsilon_symm_saddle < epsilon_2_saddle))
10213 else if (epsilon_1_saddle < epsilon_2_saddle)
10220 if (
mod(Z, 2.0) == 0)
10221 r_e_o = std::pow(10.0, -0.0170 * (E_saddle_scission + Eld) * (E_saddle_scission + Eld));
10239 if (rmode < Yasymm1)
10241 else if ((rmode > Yasymm1) && (rmode < Yasymm))
10251 N1mean = Nheavy1_eff;
10252 N1width = wNasymm1;
10258 N1mean = Nheavy2_eff;
10259 N1width = wNasymm2;
10276 while (N1r < 5.0 || N2r < 5.0)
10280 N1r =
gausshaz(0, N1mean, N1width);
10287 Z1UCD = Z /
N * N1r;
10288 Z2UCD = Z /
N * N2r;
10297 E_scission_pre =
max(epsilon_1_scission, 1.0);
10300 if (N1mean >
N * 0.50)
10314 E_scission_pre =
max(epsilon_2_scission, 1.0);
10315 if (N1mean >
N * 0.50)
10317 beta1 = (N1r - 92.0) * 0.030 + 0.60;
10319 beta1gs = ecld->beta2[
idint(N1r)][
idint(Z1UCD)];
10320 beta2gs = ecld->beta2[
idint(N2r)][
idint(Z2UCD)];
10322 beta1 =
max(beta1, beta1gs);
10323 beta2 = 1.0 - beta1;
10324 beta2 =
max(beta2, beta2gs);
10329 beta1gs = ecld->beta2[
idint(N1r)][
idint(Z1UCD)];
10330 beta2gs = ecld->beta2[
idint(N2r)][
idint(Z2UCD)];
10332 beta2 = (N2r - 92.0) * 0.030 + 0.60;
10333 beta2 =
max(beta2, beta2gs);
10334 beta1 = 1.0 - beta2;
10335 beta1 =
max(beta1, beta1gs);
10348 betags = ecld->beta2[
idint(Nsymm)][
idint(Zsymm)];
10349 beta1gs = ecld->beta2[
idint(N1r)][
idint(Z1UCD)];
10350 beta2gs = ecld->beta2[
idint(N2r)][
idint(Z2UCD)];
10351 beta =
max(0.177963 + 0.0153241 * Zsymm - 1.62037e-4 * Zsymm * Zsymm, betags);
10352 beta1 =
max(0.177963 + 0.0153241 * Z1UCD - 1.62037e-4 * Z1UCD * Z1UCD, beta1gs);
10353 beta2 =
max(0.177963 + 0.0153241 * Z2UCD - 1.62037e-4 * Z2UCD * Z2UCD, beta2gs);
10355 E_asym =
frldm(Z1UCD, N1r, beta1) +
frldm(Z2UCD, N2r, beta2) +
10356 ecoul(Z1UCD, N1r, beta1, Z2UCD, N2r, beta2, 2.0) - 2.0 *
frldm(Zsymm, Nsymm, beta) -
10357 ecoul(Zsymm, Nsymm, beta, Zsymm, Nsymm, beta, 2.0);
10358 E_scission_pre =
max(epsilon_symm_scission - E_asym, 1.);
10367 if (E_scission_pre > 5. && NbLam0 < 1)
10370 A, Z, E_scission_pre, &E_scission_post, &A_scission, &Z_scission, vx_eva_sc, vy_eva_sc, vz_eva_sc, &NbLam0);
10371 N_scission = A_scission - Z_scission;
10377 E_scission_post = E_scission_pre;
10378 N_scission = A_scission - Z_scission;
10384 N1r = N1r * N_scission /
N;
10385 N2r = N2r * N_scission /
N;
10386 Z1UCD = Z1UCD * Z_scission / Z;
10387 Z2UCD = Z2UCD * Z_scission / Z;
10424 CZ = (
frldm(Z1UCD - 1.0, N1r + 1.0, beta1) +
frldm(Z2UCD + 1.0, N2r - 1.0, beta2) +
10425 frldm(Z1UCD + 1.0, N1r - 1.0, beta1) +
frldm(Z2UCD - 1.0, N2r + 1.0, beta2) +
10426 ecoul(Z1UCD - 1.0, N1r + 1.0, beta1, Z2UCD + 1.0, N2r - 1.0, beta2, 2.0) +
10427 ecoul(Z1UCD + 1.0, N1r - 1.0, beta1, Z2UCD - 1.0, N2r + 1.0, beta2, 2.0) -
10428 2.0 *
ecoul(Z1UCD, N1r, beta1, Z2UCD, N2r, beta2, 2.0) - 2.0 *
frldm(Z1UCD, N1r, beta1) -
10429 2.0 *
frldm(Z2UCD, N2r, beta2)) *
10432 if (1.0 / A_levdens * E_scission_post < 0.0)
10433 std::cout <<
"DSQRT 1 < 0" << A_levdens <<
" " << E_scission_post << std::endl;
10435 if (0.50 * std::sqrt(1.0 / A_levdens * E_scission_post) / CZ < 0.0)
10437 std::cout <<
"DSQRT 2 < 0 " << CZ << std::endl;
10438 std::cout <<
"This event was not considered" << std::endl;
10442 ZA1width = std::sqrt(0.5 * std::sqrt(1.0 / A_levdens * E_scission_post) / CZ);
10453 ZA1width =
max(ZA1width, sigZmin);
10455 if (imode == 1 && cpol1 != 0.0)
10460 Z1rr = Z1UCD - cpol1 * A_scission / N_scission;
10463 Z1r =
gausshaz(0, Z1rr, ZA1width);
10467 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN "
10468 "CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED"
10472 if ((
utilabs(Z1rr - Z1r) > 3.0 * ZA1width) || Z1r < 1.0)
10478 if (imode == 2 && cpol2 != 0.0)
10483 Z1rr = Z1UCD - cpol2 * A_scission / N_scission;
10485 Z1r =
gausshaz(0, Z1rr, ZA1width);
10489 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN "
10490 "CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED"
10494 if ((
utilabs(Z1rr - Z1r) > 3.0 * ZA1width) || Z1r < 1.0)
10506 re1 =
frldm(Z1UCD - 1.0, N1r + 1.0, beta1) +
frldm(Z2UCD + 1.0, N2r - 1.0, beta2) +
10507 ecoul(Z1UCD - 1.0, N1r + 1.0, beta1, Z2UCD + 1.0, N2r - 1.0, beta2, d);
10508 re2 =
frldm(Z1UCD, N1r, beta1) +
frldm(Z2UCD, N2r, beta2) +
10509 ecoul(Z1UCD, N1r, beta1, Z2UCD, N2r, beta2, d);
10510 re3 =
frldm(Z1UCD + 1.0, N1r - 1.0, beta1) +
frldm(Z2UCD - 1.0, N2r + 1.0, beta2) +
10511 ecoul(Z1UCD + 1.0, N1r - 1.0, beta1, Z2UCD - 1.0, N2r + 1.0, beta2, d);
10512 eps2 = (re1 - 2.0 * re2 + re3) / 2.0;
10513 eps1 = (re3 - re1) / 2.0;
10514 DN1_POL = -eps1 / (2.0 * eps2);
10516 Z1rr = Z1UCD + DN1_POL;
10523 DN1_POL = DN1_POL - 0.6 *
Uwash(E_scission_post, Ecrit, FREDSHELL, gamma);
10524 Z1rr = Z1UCD + DN1_POL;
10530 DN1_POL = DN1_POL + 0.60 *
Uwash(E_scission_post, Ecrit, FREDSHELL, gamma);
10531 Z1rr = Z1UCD + DN1_POL;
10540 Z1r =
gausshaz(0, Z1rr, ZA1width);
10544 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN "
10545 "CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED"
10550 if ((
utilabs(Z1rr - Z1r) > 3.0 * ZA1width) || (Z1r < 1.0))
10562 z2 =
dint(Z_scission) - z1;
10564 N2 =
dint(N_scission) - N1;
10568 if ((z1 < 0) || (z2 < 0) || (a1 < 0) || (a2 < 0))
10570 std::cout <<
" -------------------------------" << std::endl;
10571 std::cout <<
" Z, A, N : " << Z <<
" " <<
A <<
" " <<
N << std::endl;
10572 std::cout << z1 <<
" " << z2 <<
" " << a1 <<
" " << a2 << std::endl;
10573 std::cout << E_scission_post <<
" " << A_levdens <<
" " << CZ << std::endl;
10575 std::cout <<
" -------------------------------" << std::endl;
10585 if (N1mean >
N * 0.50)
10589 beta2gs = ecld->beta2[
idint(N2)][
idint(z2)];
10590 if (beta2 < beta2gs)
10592 E1exc = E_scission_pre * a1 /
A + E_defo;
10593 E_defo =
frldm(z2, N2, beta2) -
frldm(z2, N2, beta2gs);
10594 E2exc = E_scission_pre * a2 /
A + E_defo;
10599 beta1gs = ecld->beta2[
idint(N1)][
idint(z1)];
10600 if (beta1 < beta1gs)
10602 E_defo =
frldm(z1, N1, beta1) -
frldm(z1, N1, beta1gs);
10603 E1exc = E_scission_pre * a1 /
A + E_defo;
10605 E2exc = E_scission_pre * a2 /
A + E_defo;
10612 if (N1mean >
N * 0.5)
10615 beta1gs = ecld->beta2[
idint(N1)][
idint(z1)];
10616 if (beta1 < beta1gs)
10618 E_defo =
frldm(z1, N1, beta1) -
frldm(z1, N1, beta1gs);
10619 E1exc = E_scission_pre * a1 /
A + E_defo;
10620 beta2gs = ecld->beta2[
idint(N2)][
idint(z2)];
10621 if (beta2 < beta2gs)
10623 E_defo =
frldm(z2, N2, beta2) -
frldm(z2, N2, beta2gs);
10624 E2exc = E_scission_pre * a2 /
A + E_defo;
10629 beta2gs = ecld->beta2[
idint(N2)][
idint(z2)];
10630 if (beta2 < beta2gs)
10632 E_defo =
frldm(z2, N2, beta2) -
frldm(z2, N2, beta2gs);
10633 E2exc = E_scission_pre * a2 /
A + E_defo;
10634 beta1gs = ecld->beta2[
idint(N1)][
idint(z1)];
10635 if (beta1 < beta1gs)
10637 E_defo =
frldm(z1, N1, beta1) -
frldm(z1, N1, beta1gs);
10638 E1exc = E_scission_pre * a1 /
A + E_defo;
10645 beta1gs = ecld->beta2[
idint(N1)][
idint(z1)];
10646 if (beta1 < beta1gs)
10648 beta2gs = ecld->beta2[
idint(N2)][
idint(z2)];
10649 if (beta2 < beta2gs)
10651 E_defo1 =
frldm(z1, N1, beta1) -
frldm(z1, N1, beta1gs);
10652 E_defo2 =
frldm(z2, N2, beta2) -
frldm(z2, N2, beta2gs);
10653 E1exc = E_scission_pre * a1 /
A + E_defo1;
10654 E2exc = E_scission_pre * a2 /
A + E_defo2;
10658 TKER = (z1 * z2 * 1.440) / (R0 * std::pow(a1, 0.333330) * (1.0 + 2.0 / 3.0 * beta1) +
10659 R0 * std::pow(a2, 0.333330) * (1.0 + 2.0 / 3.0 * beta2) + 2.0);
10661 EkinR1 = TKER * a2 /
A;
10662 EkinR2 = TKER * a1 /
A;
10663 v1 = std::sqrt(EkinR1 / a1) * 1.3887;
10664 v2 = std::sqrt(EkinR2 / a2) * 1.3887;
10667 E1exc_sigma = 5.50;
10668 E2exc_sigma = 5.50;
10672 e1 =
gausshaz(0, E1exc, E1exc_sigma);
10677 e2 =
gausshaz(0, E2exc, E2exc_sigma);
10681 (*NbLam0_par) = NbLam0;
10975 G4int INMIN, INMAX, NDIF = 0, IMEM;
10976 G4int NEVA = 0, PEVA = 0;
10986 BU_TAB_TEMP[i][0] = 0.0;
10987 BU_TAB_TEMP[i][1] = 0.0;
10988 BU_TAB_TEMP[i][2] = 0.0;
10989 BU_TAB_TEMP[i][3] = 0.0;
10990 BU_TAB_TEMP[i][4] = 0.0;
10997 if (AFP == 0 && ZFP == 0)
11002 if ((AFP == 1 && ZFP == 0) || (AFP == 1 && ZFP == 1) || (AFP == 2 && ZFP == 1) || (AFP == 3 && ZFP == 1) ||
11003 (AFP == 3 && ZFP == 2) || (AFP == 4 && ZFP == 2) || (AFP == 6 && ZFP == 2) || (AFP == 8 && ZFP == 2))
11011 if ((AFP - ZFP) == 0 && ZFP > 1)
11013 for (
G4int I = 0; I <= AFP - 2; I++)
11028 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11029 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11030 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11031 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11032 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11033 *ILOOP = *ILOOP + 1;
11051 for (
G4int I = 1; I <= 10; I++)
11065 for (
G4int I = 0; I < IMEM; I++)
11069 G4double(NDIF + ZFP + IMEM - I - 1),
11080 BU_TAB_TEMP[I + 1 + *ILOOP][0] = 1.0;
11081 BU_TAB_TEMP[I + 1 + *ILOOP][1] = 1.0;
11082 BU_TAB_TEMP[I + 1 + *ILOOP][2] = VP2X;
11083 BU_TAB_TEMP[I + 1 + *ILOOP][3] = VP2Y;
11084 BU_TAB_TEMP[I + 1 + *ILOOP][4] = VP2Z;
11089 *ILOOP = *ILOOP + IMEM;
11094 NEVA = NDIF - INMAX;
11097 for (
G4int I = 0; I < NEVA; I++)
11113 BU_TAB_TEMP[*ILOOP][0] = 0.0;
11114 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11115 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11116 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11117 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11118 *ILOOP = *ILOOP + 1;
11125 if ((AFP >= 2) && (ZFP == 0))
11127 for (
G4int I = 0; I <= AFP - 2; I++)
11143 BU_TAB_TEMP[*ILOOP][0] = 0.0;
11144 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11145 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11146 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11147 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11148 *ILOOP = *ILOOP + 1;
11161 std::cout <<
"WARNING - BU UNSTABLE: AF < ZF" << std::endl;
11166 if ((AFP >= 4) && (ZFP == 1))
11170 for (
G4int I = 0; I < AFP - 3; I++)
11186 BU_TAB_TEMP[*ILOOP][0] = 0.0;
11187 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11188 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11189 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11190 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11191 *ILOOP = *ILOOP + 1;
11203 if ((AFP == 4) && (ZFP == 3))
11210 unstable_tke(4.0, 3.0, 3.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11212 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11213 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11214 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11215 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11216 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11217 *ILOOP = *ILOOP + 1;
11219 if ((AFP == 5) && (ZFP == 2))
11226 unstable_tke(5.0, 2.0, 4.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11227 BU_TAB_TEMP[*ILOOP][0] = 0.0;
11228 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11229 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11230 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11231 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11232 *ILOOP = *ILOOP + 1;
11235 if ((AFP == 5) && (ZFP == 3))
11242 unstable_tke(5.0, 3.0, 4.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11243 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11244 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11245 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11246 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11247 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11248 *ILOOP = *ILOOP + 1;
11251 if ((AFP == 6) && (ZFP == 4))
11259 unstable_tke(6.0, 4.0, 5.0, 3.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11260 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11261 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11262 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11263 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11264 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11265 *ILOOP = *ILOOP + 1;
11271 unstable_tke(5.0, 3.0, 4.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11272 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11273 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11274 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11275 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11276 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11277 *ILOOP = *ILOOP + 1;
11279 if ((AFP == 7) && (ZFP == 2))
11286 unstable_tke(7.0, 2.0, 6.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11287 BU_TAB_TEMP[*ILOOP][0] = 0.0;
11288 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11289 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11290 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11291 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11292 *ILOOP = *ILOOP + 1;
11295 if ((AFP == 7) && (ZFP == 5))
11298 for (
int I = 0; I <= AFP - 5; I++)
11302 double(AFP - I - 1),
11303 double(ZFP - I - 1),
11313 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11314 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11315 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11316 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11317 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11318 *ILOOP = *ILOOP + 1;
11329 if ((AFP == 8) && (ZFP == 4))
11335 unstable_tke(8.0, 4.0, 4.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11336 BU_TAB_TEMP[*ILOOP][0] = 2.0;
11337 BU_TAB_TEMP[*ILOOP][1] = 4.0;
11338 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11339 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11340 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11341 *ILOOP = *ILOOP + 1;
11343 if ((AFP == 8) && (ZFP == 6))
11351 unstable_tke(8.0, 6.0, 7.0, 5.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11352 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11353 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11354 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11355 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11356 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11357 *ILOOP = *ILOOP + 1;
11362 unstable_tke(7.0, 5.0, 6.0, 4.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11363 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11364 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11365 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11366 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11367 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11368 *ILOOP = *ILOOP + 1;
11374 if ((AFP == 9) && (ZFP == 2))
11382 unstable_tke(9.0, 2.0, 8.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11383 BU_TAB_TEMP[*ILOOP][0] = 0.0;
11384 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11385 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11386 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11387 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11388 *ILOOP = *ILOOP + 1;
11394 if ((AFP == 9) && (ZFP == 5))
11401 unstable_tke(9.0, 5.0, 8.0, 4.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11402 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11403 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11404 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11405 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11406 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11407 *ILOOP = *ILOOP + 1;
11412 unstable_tke(8.0, 4.0, 4.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11413 BU_TAB_TEMP[*ILOOP][0] = 2.0;
11414 BU_TAB_TEMP[*ILOOP][1] = 4.0;
11415 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11416 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11417 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11418 *ILOOP = *ILOOP + 1;
11424 if ((AFP == 10) && (ZFP == 2))
11432 unstable_tke(10.0, 2.0, 9.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11433 BU_TAB_TEMP[*ILOOP][0] = 0.0;
11434 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11435 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11436 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11437 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11438 *ILOOP = *ILOOP + 1;
11444 unstable_tke(9.0, 2.0, 8.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11445 BU_TAB_TEMP[*ILOOP][0] = 0.0;
11446 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11447 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11448 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11449 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11450 *ILOOP = *ILOOP + 1;
11455 if ((AFP == 10) && (ZFP == 3))
11462 unstable_tke(10.0, 3.0, 9.0, 3.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11463 BU_TAB_TEMP[*ILOOP][0] = 0.0;
11464 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11465 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11466 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11467 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11468 *ILOOP = *ILOOP + 1;
11473 if ((AFP == 10) && (ZFP == 7))
11480 unstable_tke(10.0, 7.0, 9.0, 6.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11481 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11482 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11483 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11484 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11485 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11486 *ILOOP = *ILOOP + 1;
11492 if ((AFP == 11) && (ZFP == 7))
11499 unstable_tke(11.0, 7.0, 10.0, 6.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11500 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11501 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11502 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11503 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11504 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11505 *ILOOP = *ILOOP + 1;
11510 if ((AFP == 12) && (ZFP == 8))
11518 unstable_tke(12.0, 8.0, 11.0, 7.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11519 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11520 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11521 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11522 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11523 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11524 *ILOOP = *ILOOP + 1;
11529 unstable_tke(11.0, 7.0, 10.0, 6.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11530 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11531 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11532 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11533 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11534 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11535 *ILOOP = *ILOOP + 1;
11540 if ((AFP == 15) && (ZFP == 9))
11547 unstable_tke(15.0, 9.0, 14.0, 8.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11548 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11549 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11550 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11551 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11552 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11553 *ILOOP = *ILOOP + 1;
11559 if ((AFP == 16) && (ZFP == 9))
11566 unstable_tke(16.0, 9.0, 15.0, 8.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11567 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11568 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11569 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11570 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11571 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11572 *ILOOP = *ILOOP + 1;
11578 if ((AFP == 16) && (ZFP == 10))
11585 unstable_tke(16.0, 10.0, 15.0, 9.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11586 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11587 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11588 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11589 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11590 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11591 *ILOOP = *ILOOP + 1;
11596 unstable_tke(15.0, 9.0, 14.0, 8.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11597 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11598 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11599 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11600 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11601 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11602 *ILOOP = *ILOOP + 1;
11607 if ((AFP == 18) && (ZFP == 11))
11614 unstable_tke(18.0, 11.0, 17.0, 10.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11615 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11616 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11617 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11618 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11619 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11620 *ILOOP = *ILOOP + 1;
11625 if ((AFP == 19) && (ZFP == 11))
11632 unstable_tke(19.0, 11.0, 18.0, 10.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
11633 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11634 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11635 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11636 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11637 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11638 *ILOOP = *ILOOP + 1;
11643 if (ZFP >= 4 && (AFP - ZFP) == 1)
11649 for (
G4int I = 0; I < NEVA; I++)
11664 BU_TAB_TEMP[*ILOOP][0] = 0.0;
11665 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11666 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11667 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11668 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11669 *ILOOP = *ILOOP + 1;
11674 for (
G4int I = 0; I < PEVA; I++)
11689 BU_TAB_TEMP[*ILOOP][0] = 1.0;
11690 BU_TAB_TEMP[*ILOOP][1] = 1.0;
11691 BU_TAB_TEMP[*ILOOP][2] = VP2X;
11692 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
11693 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
11694 *ILOOP = *ILOOP + 1;