99 if (nucleusS > 0) nucleusS = 0;
101 G4int NbLam0 = std::abs(nucleusS);
103 Ainit = -1 * nucleusA;
104 Zinit = -1 * nucleusZ;
105 Sinit = -1 * nucleusS;
109 G4int ZFP1 = 0, AFP1 = 0, AFPIMF = 0, ZFPIMF = 0, ZFP2 = 0, AFP2 = 0, SFP1 = 0, SFP2 = 0,
111 G4double vx_eva = 0.0, vy_eva = 0.0, vz_eva = 0.0;
112 G4double VX_PREF = 0., VY_PREF = 0., VZ_PREF = 00, VP1X, VP1Y, VP1Z, VXOUT, VYOUT, VZOUT, V_CM[3],
113 VFP1_CM[3], VFP2_CM[3], VIMF_CM[3], VX2OUT, VY2OUT, VZ2OUT;
114 G4double zf = 0.0, af = 0.0, mtota = 0.0, tkeimf = 0.0, jprf0 = 0.;
115 G4int ff = 0, afpnew = 0, zfpnew = 0, aprfp = 0, zprfp = 0, IOUNSTABLE = 0, ILOOP = 0,
116 IEV_TAB = 0, IEV_TAB_TEMP = 0;
117 G4int fimf = 0, INMIN = 0, INMAX = 0;
119 G4int inum = eventnumber;
121 opt->optimfallowed = 1;
131 opt->nblan0 = NbLam0;
145 G4double T_init = 0., T_diff = 0., a_tilda = 0., a_tilda_BU = 0., EE_diff = 0., EINCL = 0.,
146 A_FINAL = 0., Z_FINAL = 0., E_FINAL = 0.;
148 G4double A_diff = 0., ASLOPE1, ASLOPE2, A_ACC, ABU_SLOPE, ABU_SUM = 0., AMEM = 0., ZMEM = 0.,
149 EMEM = 0., JMEM = 0., PX_BU_SUM = 0.0, PY_BU_SUM = 0.0, PZ_BU_SUM = 0.0, ETOT_SUM = 0.,
150 P_BU_SUM = 0., ZBU_SUM = 0., Z_Breakup_sum = 0., A_Breakup, Z_Breakup, N_Breakup, G_SYMM,
151 CZ, Sigma_Z, Z_Breakup_Mean, ZTEMP = 0., ATEMP = 0.;
153 G4double ETOT_PRF = 0.0, PXPRFP = 0., PYPRFP = 0., PZPRFP = 0., PPRFP = 0., VX1_BU = 0.,
154 VY1_BU = 0., VZ1_BU = 0., VBU2 = 0., GAMMA_REL = 1.0, Eexc_BU_SUM = 0., VX_BU_SUM = 0.,
155 VY_BU_SUM = 0., VZ_BU_SUM = 0., E_tot_BU = 0., EKIN_BU = 0., ZIMFBU = 0., AIMFBU = 0.,
156 ZFFBU = 0., AFFBU = 0., AFBU = 0., ZFBU = 0., EEBU = 0., TKEIMFBU = 0., vx_evabu = 0.,
157 vy_evabu = 0., vz_evabu = 0., Bvalue_BU = 0., P_BU = 0., ETOT_BU = 1., PX_BU = 0.,
158 PY_BU = 0., PZ_BU = 0., VX2_BU = 0., VY2_BU = 0., VZ2_BU = 0.;
160 G4int ABU_DIFF, ZBU_DIFF, NBU_DIFF;
161 G4int INEWLOOP = 0, ILOOPBU = 0;
168 std::cout <<
"Error - Remnant with a mass number A below 1." << std::endl;
173 for (
G4int j = 0; j < 3; j++) {
181 for (
G4int I2 = 0; I2 < 12; I2++)
182 BU_TAB[I1][I2] = 0.0;
183 for (
G4int I2 = 0; I2 < 6; I2++) {
184 BU_TAB_TEMP[I1][I2] = 0.0;
185 BU_TAB_TEMP1[I1][I2] = 0.0;
186 EV_TAB_TEMP[I1][I2] = 0.0;
187 EV_TAB[I1][I2] = 0.0;
188 EV_TAB_SSC[I1][I2] = 0.0;
189 EV_TEMP[I1][I2] = 0.0;
216 G4double pincl = std::sqrt(pxrem * pxrem + pyrem * pyrem + pzrem * pzrem);
218 G4double ETOT_incl = std::sqrt(pincl * pincl + (AAINCL * amu) * (AAINCL * amu));
219 G4double VX_incl =
C * pxrem / ETOT_incl;
220 G4double VY_incl =
C * pyrem / ETOT_incl;
221 G4double VZ_incl =
C * pzrem / ETOT_incl;
253 (T_freeze_out_in >= 0.0) ? T_freeze_out_in :
max(9.33 * std::exp(-0.00282 * AAINCL), 5.5);
256 ald->av * aprf + ald->as * std::pow(aprf, 2.0 / 3.0) + ald->ak * std::pow(aprf, 1.0 / 3.0);
258 T_init = std::sqrt(EINCL / a_tilda);
260 T_diff = T_init - T_freeze_out;
262 if (T_diff > 0.1 && zprf > 2. && (aprf - zprf) > 0.) {
267 for (
G4int i = 0; i < 5; i++) {
268 EE_diff = EINCL - a_tilda * T_freeze_out * T_freeze_out;
275 A_diff =
dint(EE_diff / (8.0 * 5.0 / T_freeze_out));
277 if (A_diff > AAINCL) A_diff = AAINCL;
279 A_FINAL = AAINCL - A_diff;
281 a_tilda = ald->av * A_FINAL + ald->as * std::pow(A_FINAL, 2.0 / 3.0)
282 + ald->ak * std::pow(A_FINAL, 1.0 / 3.0);
283 E_FINAL = a_tilda * T_freeze_out * T_freeze_out;
286 EE_diff = EINCL - E_FINAL;
298 Z_FINAL =
dint(zprf * A_FINAL / (aprf));
300 if (E_FINAL < 0.0) E_FINAL = 0.0;
306 A_diff = AAINCL - aprf;
316 else if (A_diff > 1.0) {
323 a_tilda = ald->av * AAINCL + ald->as * std::pow(AAINCL, 2.0 / 3.0)
324 + ald->ak * std::pow(AAINCL, 1.0 / 3.0);
326 E_FINAL = a_tilda * T_freeze_out * T_freeze_out;
329 (ASLOPE1 - ASLOPE2) / 4.0 * (E_FINAL / AAINCL) + ASLOPE1 - (ASLOPE1 - ASLOPE2) * 7.0 / 4.0;
341 if (ABU_SLOPE > -1.01) ABU_SLOPE = -1.01;
344 Z_Breakup_sum = Z_FINAL;
348 for (
G4int i = 0; i < 100; i++) {
355 std::cout <<
"WARNING: IPOWERLIMHAZ CALLED MORE THAN 100 TIMES WHEN "
356 "CALCULATING A_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED: "
357 << A_Breakup << std::endl;
361 if (A_Breakup > AAINCL)
goto mult4326;
363 if (A_Breakup <= 0.0) {
364 std::cout <<
"A_BREAKUP <= 0 " << std::endl;
368 A_ACC = A_ACC + A_Breakup;
370 if (A_ACC <= A_diff) {
371 Z_Breakup_Mean =
dint(A_Breakup * ZAINCL / AAINCL);
373 Z_Breakup_sum = Z_Breakup_sum + Z_Breakup_Mean;
376 G_SYMM = 34.2281 - 5.14037 * E_FINAL / AAINCL;
377 if (E_FINAL / AAINCL < 2.0) G_SYMM = 25.0;
378 if (E_FINAL / AAINCL > 4.0) G_SYMM = 15.0;
383 CZ = 2.0 * G_SYMM * 4.0 / A_Breakup;
387 Sigma_Z = std::sqrt(T_freeze_out / CZ);
395 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN "
396 "CALCULATING Z_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE "
398 << A_Breakup <<
" " << Z_Breakup << std::endl;
402 if (Z_Breakup < 0.0)
goto mult4333;
403 if ((A_Breakup - Z_Breakup) < 0.0)
goto mult4333;
404 if ((A_Breakup - Z_Breakup) == 0.0 && Z_Breakup != 1.0)
goto mult4333;
406 if (Z_Breakup >= ZAINCL) {
409 std::cout <<
"Z_BREAKUP RESAMPLED MORE THAN 10 TIMES; EVENT WILL "
410 "BE RESAMPLED AGAIN "
420 if (Z_Breakup > 2.0) {
421 if (
idnint(A_Breakup - Z_Breakup) < INMIN
422 ||
idnint(A_Breakup - Z_Breakup) > (INMAX + 5))
435 N_Breakup = A_Breakup - Z_Breakup;
436 BU_TAB[I_Breakup][0] =
dint(Z_Breakup);
437 BU_TAB[I_Breakup][1] =
dint(A_Breakup);
438 ABU_SUM = ABU_SUM + BU_TAB[i][1];
439 ZBU_SUM = ZBU_SUM + BU_TAB[i][0];
442 BU_TAB[I_Breakup][3] = 0.0;
443 I_Breakup = I_Breakup + 1;
444 IMULTBU = IMULTBU + 1;
463 ABU_DIFF =
idnint(ABU_SUM + aprf - AAINCL);
464 ZBU_DIFF =
idnint(ZBU_SUM + zprf - ZAINCL);
465 NBU_DIFF =
idnint((ABU_SUM - ZBU_SUM) + (aprf - zprf) - (AAINCL - ZAINCL));
467 if (IMULTBU > 200) std::cout <<
"WARNING - MORE THAN 200 BU " << IMULTBU << std::endl;
469 if (IMULTBU < 1) std::cout <<
"WARNING - LESS THAN 1 BU " << IMULTBU << std::endl;
473 for (
G4int i = 0; i < IMULTBU; i++)
476 while (NBU_DIFF != 0 && ZBU_DIFF != 0) {
486 std::cout <<
"WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING "
487 "N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED."
492 if (IMEM_BU[IEL] == 1)
goto mult5555;
493 if (!(IEL < 200)) std::cout <<
"5555:" << IEL << RHAZ << IMULTBU << std::endl;
494 if (IEL < 0) std::cout <<
"5555:" << IEL << RHAZ << IMULTBU << std::endl;
495 if (IEL <= IMULTBU) {
496 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0] -
DSIGN(1.0,
G4double(NBU_DIFF)));
498 else if (IEL > IMULTBU) {
501 if (N_Breakup < 0.0) {
505 if (IEL <= IMULTBU) {
508 else if (IEL > IMULTBU) {
515 if (ZTEMP < 1.0 && N_Breakup < 1.0) {
526 if (IEL <= IMULTBU) {
527 BU_TAB[IEL][0] =
dint(ZTEMP);
528 BU_TAB[IEL][1] =
dint(ZTEMP + N_Breakup);
530 else if (IEL > IMULTBU) {
532 aprf =
dint(ZTEMP + N_Breakup);
534 NBU_DIFF = NBU_DIFF -
ISIGN(1, NBU_DIFF);
535 ZBU_DIFF = ZBU_DIFF -
ISIGN(1, ZBU_DIFF);
539 for (
G4int i = 0; i < IMULTBU; i++)
542 if (NBU_DIFF != 0 && ZBU_DIFF == 0) {
543 while (NBU_DIFF > 0 || NBU_DIFF < 0) {
549 std::cout <<
"WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING "
550 "N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED."
555 if (IMEM_BU[IEL] == 1)
goto mult5556;
557 if (IPROBA > IMULTBU + 1 && NBU_DIFF > 0) {
558 std::cout <<
"###',IPROBA,IMULTBU,NBU_DIFF,ZBU_DIFF,T_freeze_out" << std::endl;
560 if (IEL <= IMULTBU) {
561 BU_TAB[IEL][1] =
dint(BU_TAB[IEL][1] -
G4double(NBU_DIFF));
564 if (IEL > IMULTBU) aprf =
dint(aprf -
G4double(NBU_DIFF));
568 if (!(IEL < 200)) std::cout <<
"5556:" << IEL << RHAZ << IMULTBU << std::endl;
569 if (IEL < 0) std::cout <<
"5556:" << IEL << RHAZ << IMULTBU << std::endl;
570 if (IEL <= IMULTBU) {
571 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0] -
DSIGN(1.0,
G4double(NBU_DIFF)));
573 else if (IEL > IMULTBU) {
576 if (N_Breakup < 0.0) {
580 if (IEL <= IMULTBU) {
581 ATEMP =
dint(BU_TAB[IEL][0] + N_Breakup);
583 else if (IEL > IMULTBU) {
584 ATEMP =
dint(zprf + N_Breakup);
586 if ((ATEMP - N_Breakup) < 1.0 && N_Breakup < 1.0) {
596 BU_TAB[IEL][1] =
dint(BU_TAB[IEL][0] + N_Breakup);
597 else if (IEL > IMULTBU)
598 aprf =
dint(zprf + N_Breakup);
600 NBU_DIFF = NBU_DIFF -
ISIGN(1, NBU_DIFF);
604 for (
G4int i = 0; i < IMULTBU; i++)
608 if (ZBU_DIFF != 0 && NBU_DIFF == 0) {
609 while (ZBU_DIFF > 0 || ZBU_DIFF < 0) {
615 std::cout <<
"WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING "
616 "N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED."
621 if (IMEM_BU[IEL] == 1)
goto mult5557;
623 if (IPROBA > IMULTBU + 1 && ZBU_DIFF > 0) {
624 std::cout <<
"###',IPROBA,IMULTBU,NBU_DIFF,ZBU_DIFF,T_freeze_out" << std::endl;
626 if (IEL <= IMULTBU) {
627 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0]);
628 BU_TAB[IEL][0] =
dint(BU_TAB[IEL][0] -
G4double(ZBU_DIFF));
629 BU_TAB[IEL][1] =
dint(BU_TAB[IEL][0] + N_Breakup);
633 N_Breakup = aprf - zprf;
635 aprf =
dint(zprf + N_Breakup);
640 if (!(IEL < 200)) std::cout <<
"5557:" << IEL << RHAZ << IMULTBU << std::endl;
641 if (IEL < 0) std::cout <<
"5557:" << IEL << RHAZ << IMULTBU << std::endl;
642 if (IEL <= IMULTBU) {
643 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0]);
646 else if (IEL > IMULTBU) {
647 N_Breakup =
dint(aprf - zprf);
650 ATEMP =
dint(ZTEMP + N_Breakup);
655 if ((ATEMP - ZTEMP) < 0.0) {
659 if ((ATEMP - ZTEMP) < 1.0 && ZTEMP < 1.0) {
663 if (IEL <= IMULTBU) {
664 BU_TAB[IEL][0] =
dint(ZTEMP);
665 BU_TAB[IEL][1] =
dint(ZTEMP + N_Breakup);
670 aprf =
dint(ZTEMP + N_Breakup);
673 ZBU_DIFF = ZBU_DIFF -
ISIGN(1, ZBU_DIFF);
683 for (
G4int i = 0; i < IMULTBU; i++) {
689 if (BU_TAB[i][0] > 2.0) {
690 a_tilda_BU = ald->av * BU_TAB[i][1] + ald->as * std::pow(BU_TAB[i][1], 2.0 / 3.0)
691 + ald->ak * std::pow(BU_TAB[i][1], 1.0 / 3.0);
692 BU_TAB[i][2] = a_tilda_BU * T_freeze_out * T_freeze_out;
698 if (BU_TAB[i][0] > ZMEM) {
708 BU_TAB[IMEM][0] = zprf;
709 BU_TAB[IMEM][1] = aprf;
710 BU_TAB[IMEM][2] = ee;
711 BU_TAB[IMEM][3] = jprf;
723 for (
G4int i = 0; i < IMULTBU; i++) {
724 ABU_SUM = ABU_SUM + BU_TAB[i][1];
725 ZBU_SUM = ZBU_SUM + BU_TAB[i][0];
727 ABU_DIFF =
idnint(ABU_SUM - AAINCL);
728 ZBU_DIFF =
idnint(ZBU_SUM - ZAINCL);
730 if (ABU_DIFF != 0 || ZBU_DIFF != 0)
731 std::cout <<
"Problem of mass in BU " << ABU_DIFF <<
" " << ZBU_DIFF << std::endl;
738 AMOMENT(AAINCL, aprf, 1, &PXPRFP, &PYPRFP, &PZPRFP);
739 PPRFP = std::sqrt(PXPRFP * PXPRFP + PYPRFP * PYPRFP + PZPRFP * PZPRFP);
742 ETOT_PRF = std::sqrt(PPRFP * PPRFP + (aprf * amu) * (aprf * amu));
743 VX_PREF =
C * PXPRFP / ETOT_PRF;
744 VY_PREF =
C * PYPRFP / ETOT_PRF;
745 VZ_PREF =
C * PZPRFP / ETOT_PRF;
748 tke_bu(zprf, aprf, ZAINCL, AAINCL, &VX1_BU, &VY1_BU, &VZ1_BU);
755 lorentz_boost(VX1_BU, VY1_BU, VZ1_BU, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT, &VZOUT);
762 VBU2 = VX_PREF * VX_PREF + VY_PREF * VY_PREF + VZ_PREF * VZ_PREF;
763 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C *
C));
764 ETOT_PRF = aprf * amu / GAMMA_REL;
765 PXPRFP = ETOT_PRF * VX_PREF /
C;
766 PYPRFP = ETOT_PRF * VY_PREF /
C;
767 PZPRFP = ETOT_PRF * VZ_PREF /
C;
781 for (I_Breakup = 0; I_Breakup < IMULTBU; I_Breakup++) {
785 Eexc_BU_SUM = Eexc_BU_SUM + BU_TAB[I_Breakup][2];
787 AMOMENT(AAINCL, BU_TAB[I_Breakup][1], 1, &PX_BU, &PY_BU, &PZ_BU);
788 P_BU = std::sqrt(PX_BU * PX_BU + PY_BU * PY_BU + PZ_BU * PZ_BU);
792 std::sqrt(P_BU * P_BU + (BU_TAB[I_Breakup][1] * amu) * (BU_TAB[I_Breakup][1] * amu));
793 BU_TAB[I_Breakup][4] =
C * PX_BU / ETOT_BU;
794 BU_TAB[I_Breakup][5] =
C * PY_BU / ETOT_BU;
795 BU_TAB[I_Breakup][6] =
C * PZ_BU / ETOT_BU;
797 tke_bu(BU_TAB[I_Breakup][0], BU_TAB[I_Breakup][1], ZAINCL, AAINCL, &VX2_BU, &VY2_BU, &VZ2_BU);
804 lorentz_boost(VX2_BU, VY2_BU, VZ2_BU, BU_TAB[I_Breakup][4], BU_TAB[I_Breakup][5],
805 BU_TAB[I_Breakup][6], &VXOUT, &VYOUT, &VZOUT);
807 BU_TAB[I_Breakup][4] = VXOUT;
808 BU_TAB[I_Breakup][5] = VYOUT;
809 BU_TAB[I_Breakup][6] = VZOUT;
812 VBU2 = BU_TAB[I_Breakup][4] * BU_TAB[I_Breakup][4]
813 + BU_TAB[I_Breakup][5] * BU_TAB[I_Breakup][5]
814 + BU_TAB[I_Breakup][6] * BU_TAB[I_Breakup][6];
815 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C *
C));
816 ETOT_BU = BU_TAB[I_Breakup][1] * amu / GAMMA_REL;
817 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
818 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
819 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
821 PX_BU_SUM = PX_BU_SUM + PX_BU;
822 PY_BU_SUM = PY_BU_SUM + PY_BU;
823 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
828 P_BU_SUM = std::sqrt(PX_BU_SUM * PX_BU_SUM + PY_BU_SUM * PY_BU_SUM + PZ_BU_SUM * PZ_BU_SUM);
831 ETOT_SUM = std::sqrt(P_BU_SUM * P_BU_SUM + (AAINCL * amu) * (AAINCL * amu));
833 VX_BU_SUM =
C * PX_BU_SUM / ETOT_SUM;
834 VY_BU_SUM =
C * PY_BU_SUM / ETOT_SUM;
835 VZ_BU_SUM =
C * PZ_BU_SUM / ETOT_SUM;
842 lorentz_boost(-VX_BU_SUM, -VY_BU_SUM, -VZ_BU_SUM, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT,
849 VBU2 = VX_PREF * VX_PREF + VY_PREF * VY_PREF + VZ_PREF * VZ_PREF;
850 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C *
C));
851 ETOT_PRF = aprf * amu / GAMMA_REL;
852 PXPRFP = ETOT_PRF * VX_PREF /
C;
853 PYPRFP = ETOT_PRF * VY_PREF /
C;
854 PZPRFP = ETOT_PRF * VZ_PREF /
C;
865 EKIN_BU = aprf * amu / GAMMA_REL - aprf * amu;
867 for (I_Breakup = 0; I_Breakup < IMULTBU; I_Breakup++) {
873 lorentz_boost(-VX_BU_SUM, -VY_BU_SUM, -VZ_BU_SUM, BU_TAB[I_Breakup][4], BU_TAB[I_Breakup][5],
874 BU_TAB[I_Breakup][6], &VXOUT, &VYOUT, &VZOUT);
876 BU_TAB[I_Breakup][4] = VXOUT;
877 BU_TAB[I_Breakup][5] = VYOUT;
878 BU_TAB[I_Breakup][6] = VZOUT;
880 VBU2 = BU_TAB[I_Breakup][4] * BU_TAB[I_Breakup][4]
881 + BU_TAB[I_Breakup][5] * BU_TAB[I_Breakup][5]
882 + BU_TAB[I_Breakup][6] * BU_TAB[I_Breakup][6];
883 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C *
C));
885 ETOT_BU = BU_TAB[I_Breakup][1] * amu / GAMMA_REL;
887 EKIN_BU = EKIN_BU + BU_TAB[I_Breakup][1] * amu / GAMMA_REL - BU_TAB[I_Breakup][1] * amu;
889 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
890 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
891 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
892 E_tot_BU = E_tot_BU + ETOT_BU;
894 PX_BU_SUM = PX_BU_SUM + PX_BU;
895 PY_BU_SUM = PY_BU_SUM + PY_BU;
896 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
899 if (std::abs(PX_BU_SUM) > 10. || std::abs(PY_BU_SUM) > 10. || std::abs(PZ_BU_SUM) > 10.) {
901 P_BU_SUM = std::sqrt(PX_BU_SUM * PX_BU_SUM + PY_BU_SUM * PY_BU_SUM + PZ_BU_SUM * PZ_BU_SUM);
904 ETOT_SUM = std::sqrt(P_BU_SUM * P_BU_SUM + (AAINCL * amu) * (AAINCL * amu));
906 VX_BU_SUM =
C * PX_BU_SUM / ETOT_SUM;
907 VY_BU_SUM =
C * PY_BU_SUM / ETOT_SUM;
908 VZ_BU_SUM =
C * PZ_BU_SUM / ETOT_SUM;
915 lorentz_boost(-VX_BU_SUM, -VY_BU_SUM, -VZ_BU_SUM, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT,
922 VBU2 = VX_PREF * VX_PREF + VY_PREF * VY_PREF + VZ_PREF * VZ_PREF;
923 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C *
C));
924 ETOT_PRF = aprf * amu / GAMMA_REL;
925 PXPRFP = ETOT_PRF * VX_PREF /
C;
926 PYPRFP = ETOT_PRF * VY_PREF /
C;
927 PZPRFP = ETOT_PRF * VZ_PREF /
C;
938 EKIN_BU = aprf * amu / GAMMA_REL - aprf * amu;
940 for (I_Breakup = 0; I_Breakup < IMULTBU; I_Breakup++) {
946 lorentz_boost(-VX_BU_SUM, -VY_BU_SUM, -VZ_BU_SUM, BU_TAB[I_Breakup][4],
947 BU_TAB[I_Breakup][5], BU_TAB[I_Breakup][6], &VXOUT, &VYOUT, &VZOUT);
949 BU_TAB[I_Breakup][4] = VXOUT;
950 BU_TAB[I_Breakup][5] = VYOUT;
951 BU_TAB[I_Breakup][6] = VZOUT;
953 VBU2 = BU_TAB[I_Breakup][4] * BU_TAB[I_Breakup][4]
954 + BU_TAB[I_Breakup][5] * BU_TAB[I_Breakup][5]
955 + BU_TAB[I_Breakup][6] * BU_TAB[I_Breakup][6];
956 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C *
C));
958 ETOT_BU = BU_TAB[I_Breakup][1] * amu / GAMMA_REL;
960 EKIN_BU = EKIN_BU + BU_TAB[I_Breakup][1] * amu / GAMMA_REL - BU_TAB[I_Breakup][1] * amu;
962 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
963 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
964 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
965 E_tot_BU = E_tot_BU + ETOT_BU;
967 PX_BU_SUM = PX_BU_SUM + PX_BU;
968 PY_BU_SUM = PY_BU_SUM + PY_BU;
969 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
978 for (
G4int i = 0; i < IMULTBU; i++) {
979 if (BU_TAB[i][0] < 3.0 || BU_TAB[i][0] == BU_TAB[i][1]) {
981 BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VP1X, &VP1Y, &VP1Z, BU_TAB_TEMP,
984 if (IOUNSTABLE > 0) {
993 for (
int IJ = 0; IJ < ILOOP; IJ++) {
994 BU_TAB[IMULTBU + INEWLOOP + IJ][0] = BU_TAB_TEMP[IJ][0];
995 BU_TAB[IMULTBU + INEWLOOP + IJ][1] = BU_TAB_TEMP[IJ][1];
996 BU_TAB[IMULTBU + INEWLOOP + IJ][4] = BU_TAB_TEMP[IJ][2];
997 BU_TAB[IMULTBU + INEWLOOP + IJ][5] = BU_TAB_TEMP[IJ][3];
998 BU_TAB[IMULTBU + INEWLOOP + IJ][6] = BU_TAB_TEMP[IJ][4];
999 BU_TAB[IMULTBU + INEWLOOP + IJ][2] = 0.0;
1000 BU_TAB[IMULTBU + INEWLOOP + IJ][3] = 0.0;
1003 INEWLOOP = INEWLOOP + ILOOP;
1010 IMULTBU = IMULTBU + INEWLOOP;
1012 opt->optimfallowed = 1;
1018 std::vector<G4double> problamb(IMULTBU);
1021 for (
G4int i = 0; i < IMULTBU; i++)
1022 sumN = sumN + BU_TAB[i][1] - BU_TAB[i][0];
1024 for (
G4int i = 0; i < IMULTBU; i++) {
1025 problamb[i] = (BU_TAB[i][1] - BU_TAB[i][0]) / sumN;
1028 std::vector<G4int> Nblamb(IMULTBU, 0);
1029 for (
G4int j = 0; j < NbLam0;) {
1030 G4double probtotal = (aprf - zprf) / sumN;
1033 if (ran <= probtotal) {
1037 for (
G4int i = 0; i < IMULTBU; i++) {
1039 if (probtotal < ran && ran <= probtotal + problamb[i]) {
1040 Nblamb[i] = Nblamb[i] + 1;
1043 probtotal = probtotal + problamb[i];
1049 for (
G4int i = 0; i < IMULTBU; i++) {
1050 EEBU = BU_TAB[i][2];
1051 BU_TAB[i][10] = BU_TAB[i][6];
1053 if (BU_TAB[i][0] > 2.0) {
1054 G4int nbl = Nblamb[i];
1055 evapora(BU_TAB[i][0], BU_TAB[i][1], &EEBU, 0.0, &ZFBU, &AFBU, &mtota, &vz_evabu, &vx_evabu,
1056 &vy_evabu, &ff, &fimf, &ZIMFBU, &AIMFBU, &TKEIMFBU, &jprfbu, &inttype, &inum,
1057 EV_TEMP, &IEV_TAB_TEMP, &nbl);
1060 BU_TAB[i][9] = jprfbu;
1064 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1065 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1066 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1067 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1073 lorentz_boost(BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], EV_TEMP[IJ][2], EV_TEMP[IJ][3],
1074 EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1075 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1076 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1077 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1079 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1087 lorentz_boost(vx_evabu, vy_evabu, vz_evabu, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6],
1088 &VXOUT, &VYOUT, &VZOUT);
1089 BU_TAB[i][4] = VXOUT;
1090 BU_TAB[i][5] = VYOUT;
1091 BU_TAB[i][6] = VZOUT;
1094 BU_TAB[i][7] =
dint(ZFBU);
1095 BU_TAB[i][8] =
dint(AFBU);
1096 BU_TAB[i][11] = nbl;
1106 opt->optimfallowed = 0;
1109 G4double EkinR1 = TKEIMFBU * AIMFBU / (AFBU + AIMFBU);
1110 G4double EkinR2 = TKEIMFBU * AFBU / (AFBU + AIMFBU);
1111 G4double V1 = std::sqrt(EkinR1 / AFBU) * 1.3887;
1112 G4double V2 = std::sqrt(EkinR2 / AIMFBU) * 1.3887;
1114 G4double VPERP1 = std::sqrt(
V1 *
V1 - VZ1_IMF * VZ1_IMF);
1116 G4double VX1_IMF = VPERP1 * std::sin(ALPHA1);
1117 G4double VY1_IMF = VPERP1 * std::cos(ALPHA1);
1122 G4double EEIMFP = EEBU * AFBU / (AFBU + AIMFBU);
1123 G4double EEIMF = EEBU * AIMFBU / (AFBU + AIMFBU);
1126 G4double IINERTTOT = 0.40 * 931.490 * 1.160 * 1.160
1127 * (std::pow(AIMFBU, 5.0 / 3.0) + std::pow(AFBU, 5.0 / 3.0))
1128 + 931.490 * 1.160 * 1.160 * AIMFBU * AFBU / (AIMFBU + AFBU)
1129 * (std::pow(AIMFBU, 1. / 3.) + std::pow(AFBU, 1. / 3.))
1130 * (std::pow(AIMFBU, 1. / 3.) + std::pow(AFBU, 1. / 3.));
1133 BU_TAB[i][9] * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(AFBU, 5.0 / 3.0) / IINERTTOT;
1135 BU_TAB[i][9] * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(AIMFBU, 5.0 / 3.0) / IINERTTOT;
1142 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VXOUT,
1144 BU_TAB[i][4] = VXOUT;
1145 BU_TAB[i][5] = VYOUT;
1146 BU_TAB[i][6] = VZOUT;
1148 G4double vx1ev_imf = 0., vy1ev_imf = 0., vz1ev_imf = 0., zdummy = 0., adummy = 0.,
1149 tkedummy = 0., jprf1 = 0.;
1154 G4double pbH = (AFBU - ZFBU) / (AFBU - ZFBU + AIMFBU - ZIMFBU);
1155 for (
G4int j = 0; j < nbl; j++) {
1164 evapora(ZFBU, AFBU, &EEIMFP, JPRFHEAVY, &ZFFBU, &AFFBU, &mtota, &vz1ev_imf, &vx1ev_imf,
1165 &vy1ev_imf, &FFBU1, &FIMFBU1, &zdummy, &adummy, &tkedummy, &jprf1, &inttype,
1166 &inum, EV_TEMP, &IEV_TAB_TEMP, &NbLamH);
1168 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1169 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1170 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1171 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1177 lorentz_boost(BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], EV_TEMP[IJ][2], EV_TEMP[IJ][3],
1178 EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1179 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1180 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1181 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1183 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1185 BU_TAB[i][7] =
dint(ZFFBU);
1186 BU_TAB[i][8] =
dint(AFFBU);
1187 BU_TAB[i][11] = NbLamH;
1192 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6],
1193 &VXOUT, &VYOUT, &VZOUT);
1194 BU_TAB[i][4] = VXOUT;
1195 BU_TAB[i][5] = VYOUT;
1196 BU_TAB[i][6] = VZOUT;
1200 opt->optimfallowed = 0;
1203 G4double zffimf, affimf, zdummy1, adummy1, tkedummy1, jprf2, vx2ev_imf, vy2ev_imf,
1206 evapora(ZIMFBU, AIMFBU, &EEIMF, JPRFLIGHT, &zffimf, &affimf, &mtota, &vz2ev_imf,
1207 &vx2ev_imf, &vy2ev_imf, &FFBU2, &FIMFBU2, &zdummy1, &adummy1, &tkedummy1, &jprf2,
1208 &inttype, &inum, EV_TEMP, &IEV_TAB_TEMP, &NbLamimf);
1210 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1211 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1212 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1213 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1220 lorentz_boost(BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], EV_TEMP[IJ][2], EV_TEMP[IJ][3],
1221 EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1222 lorentz_boost(VX2_IMF, VY2_IMF, VZ2_IMF, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1224 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
1225 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
1226 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
1228 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1230 BU_TAB[IMULTBU + ILOOPBU][0] = BU_TAB[i][0];
1231 BU_TAB[IMULTBU + ILOOPBU][1] = BU_TAB[i][1];
1232 BU_TAB[IMULTBU + ILOOPBU][2] = BU_TAB[i][2];
1233 BU_TAB[IMULTBU + ILOOPBU][3] = BU_TAB[i][3];
1234 BU_TAB[IMULTBU + ILOOPBU][7] =
dint(zffimf);
1235 BU_TAB[IMULTBU + ILOOPBU][8] =
dint(affimf);
1236 BU_TAB[IMULTBU + ILOOPBU][11] = NbLamimf;
1238 lorentz_boost(VX2_IMF, VY2_IMF, VZ2_IMF, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VXOUT,
1240 lorentz_boost(vx2ev_imf, vy2ev_imf, vz2ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1242 BU_TAB[IMULTBU + ILOOPBU][4] = VX2OUT;
1243 BU_TAB[IMULTBU + ILOOPBU][5] = VY2OUT;
1244 BU_TAB[IMULTBU + ILOOPBU][6] = VZ2OUT;
1245 ILOOPBU = ILOOPBU + 1;
1253 BU_TAB[i][7] = BU_TAB[i][0];
1254 BU_TAB[i][8] = BU_TAB[i][1];
1258 BU_TAB[i][11] = Nblamb[i];
1262 IMULTBU = IMULTBU + ILOOPBU;
1270 for (
G4int i = 0; i < IMULTBU; i++) {
1271 ABU_SUM = ABU_SUM + BU_TAB[i][8];
1272 ZBU_SUM = ZBU_SUM + BU_TAB[i][7];
1274 BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VP1X, &VP1Y, &VP1Z, BU_TAB_TEMP1,
1282 if (IOUNSTABLE > 0) {
1284 ABU_SUM = ABU_SUM +
G4double(afpnew) - BU_TAB[i][8];
1285 ZBU_SUM = ZBU_SUM +
G4double(zfpnew) - BU_TAB[i][7];
1288 BU_TAB[i][4] = VP1X;
1289 BU_TAB[i][5] = VP1Y;
1290 BU_TAB[i][6] = VP1Z;
1293 for (
G4int IJ = 0; IJ < ILOOP; IJ++) {
1294 BU_TAB[IMULTBU + INEWLOOP + IJ][7] = BU_TAB_TEMP1[IJ][0];
1295 BU_TAB[IMULTBU + INEWLOOP + IJ][8] = BU_TAB_TEMP1[IJ][1];
1296 BU_TAB[IMULTBU + INEWLOOP + IJ][4] = BU_TAB_TEMP1[IJ][2];
1297 BU_TAB[IMULTBU + INEWLOOP + IJ][5] = BU_TAB_TEMP1[IJ][3];
1298 BU_TAB[IMULTBU + INEWLOOP + IJ][6] = BU_TAB_TEMP1[IJ][4];
1299 BU_TAB[IMULTBU + INEWLOOP + IJ][2] = 0.0;
1300 BU_TAB[IMULTBU + INEWLOOP + IJ][3] = 0.0;
1301 BU_TAB[IMULTBU + INEWLOOP + IJ][0] = BU_TAB[i][0];
1302 BU_TAB[IMULTBU + INEWLOOP + IJ][1] = BU_TAB[i][1];
1303 BU_TAB[IMULTBU + INEWLOOP + IJ][11] = BU_TAB[i][11];
1304 ABU_SUM = ABU_SUM + BU_TAB[IMULTBU + INEWLOOP + IJ][8];
1305 ZBU_SUM = ZBU_SUM + BU_TAB[IMULTBU + INEWLOOP + IJ][7];
1308 INEWLOOP = INEWLOOP + ILOOP;
1313 IMULTBU = IMULTBU + INEWLOOP;
1316 lorentz_boost(VX_incl, VY_incl, VZ_incl, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT, &VZOUT);
1321 for (
G4int i = 0; i < IMULTBU; i++) {
1322 lorentz_boost(VX_incl, VY_incl, VZ_incl, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VXOUT,
1324 BU_TAB[i][4] = VXOUT;
1325 BU_TAB[i][5] = VYOUT;
1326 BU_TAB[i][6] = VZOUT;
1328 for (
G4int i = 0; i < IEV_TAB; i++) {
1329 lorentz_boost(VX_incl, VY_incl, VZ_incl, EV_TAB[i][2], EV_TAB[i][3], EV_TAB[i][4], &VXOUT,
1331 EV_TAB[i][2] = VXOUT;
1332 EV_TAB[i][3] = VYOUT;
1333 EV_TAB[i][4] = VZOUT;
1335 if (IMULTBU > 200) std::cout <<
"IMULTBU>200 " << IMULTBU << std::endl;
1344 if (IMULTIFR == 0) {
1352 if (IMULTIFR == 1) {
1358 opt->optimfallowed = 1;
1365 if (zprfp <= 2 && zprfp < aprfp) {
1383 if (zprfp <= 2 && zprfp == aprfp) {
1384 unstable_nuclei(aprfp, zprfp, &afpnew, &zfpnew, IOUNSTABLE, VX_PREF, VY_PREF, VZ_PREF, &VP1X,
1385 &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
1391 for (
G4int I = 0; I < ILOOP; I++) {
1392 for (
G4int IJ = 0; IJ < 6; IJ++)
1393 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1395 IEV_TAB = IEV_TAB + ILOOP;
1411 if (zprfp == aprfp) {
1412 unstable_nuclei(aprfp, zprfp, &afpnew, &zfpnew, IOUNSTABLE, VX_PREF, VY_PREF, VZ_PREF, &VP1X,
1413 &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
1419 for (
G4int I = 0; I < ILOOP; I++) {
1420 for (
G4int IJ = 0; IJ < 6; IJ++)
1421 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1423 IEV_TAB = IEV_TAB + ILOOP;
1438 evapora(zprf, aprf, &ee, jprf, &zf, &af, &mtota, &vz_eva, &vx_eva, &vy_eva, &ff, &fimf, &zimf,
1439 &aimf, &tkeimf, &jprf0, &inttype, &inum, EV_TEMP, &IEV_TAB_TEMP, &NbLam0);
1441 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1442 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1443 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1444 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1450 lorentz_boost(VX_PREF, VY_PREF, VZ_PREF, EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4], &VXOUT,
1452 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1453 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1454 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1456 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1461 lorentz_boost(VX_PREF, VY_PREF, VZ_PREF, vx_eva, vy_eva, vz_eva, &VXOUT, &VYOUT, &VZOUT);
1466 if (ff == 0 && fimf == 0) {
1478 VFP1_CM[0] = V_CM[0];
1479 VFP1_CM[1] = V_CM[1];
1480 VFP1_CM[2] = V_CM[2];
1481 for (
G4int j = 0; j < 3; j++) {
1487 if (ff == 1 && fimf == 0) ftype = 1;
1488 if (ff == 0 && fimf == 1)
1499 if (NbLam0 > 0) varntp->kfis = 20;
1502 G4int IEV_TAB_FIS = 0, imode = 0;
1504 G4double vx1_fission = 0., vy1_fission = 0., vz1_fission = 0.;
1505 G4double vx2_fission = 0., vy2_fission = 0., vz2_fission = 0.;
1506 G4double vx_eva_sc = 0., vy_eva_sc = 0., vz_eva_sc = 0.;
1508 fission(af, zf, ee, jprf0, &vx1_fission, &vy1_fission, &vz1_fission, &vx2_fission, &vy2_fission,
1509 &vz2_fission, &ZFP1, &AFP1, &SFP1, &ZFP2, &AFP2, &SFP2, &imode, &vx_eva_sc, &vy_eva_sc,
1510 &vz_eva_sc, EV_TEMP, &IEV_TAB_FIS, &NbLam0);
1512 for (
G4int IJ = 0; IJ < IEV_TAB_FIS; IJ++) {
1513 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1514 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1515 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1521 lorentz_boost(V_CM[0], V_CM[1], V_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4],
1522 &VXOUT, &VYOUT, &VZOUT);
1523 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1524 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1525 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1527 IEV_TAB = IEV_TAB + IEV_TAB_FIS;
1543 lorentz_boost(vx1_fission, vy1_fission, vz1_fission, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT,
1545 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1546 VFP1_CM[0] = VX2OUT;
1547 VFP1_CM[1] = VY2OUT;
1548 VFP1_CM[2] = VZ2OUT;
1555 lorentz_boost(vx2_fission, vy2_fission, vz2_fission, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT,
1557 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1558 VFP2_CM[0] = VX2OUT;
1559 VFP2_CM[1] = VY2OUT;
1560 VFP2_CM[2] = VZ2OUT;
1566 else if (ftype == 2) {
1571 opt->optimfallowed = 1;
1576 G4double pbH = (af - zf) / (af - zf + aimf - zimf);
1578 for (
G4int i = 0; i < NbLam0; i++) {
1588 G4double EkinR1 = tkeimf * aimf / (af + aimf);
1589 G4double EkinR2 = tkeimf * af / (af + aimf);
1590 G4double V1 = std::sqrt(EkinR1 / af) * 1.3887;
1591 G4double V2 = std::sqrt(EkinR2 / aimf) * 1.3887;
1593 G4double VPERP1 = std::sqrt(
V1 *
V1 - VZ1_IMF * VZ1_IMF);
1595 G4double VX1_IMF = VPERP1 * std::sin(ALPHA1);
1596 G4double VY1_IMF = VPERP1 * std::cos(ALPHA1);
1601 G4double EEIMFP = ee * af / (af + aimf);
1602 G4double EEIMF = ee * aimf / (af + aimf);
1606 0.40 * 931.490 * 1.160 * 1.160 * (std::pow(aimf, 5.0 / 3.0) + std::pow(af, 5.0 / 3.0))
1607 + 931.490 * 1.160 * 1.160 * aimf * af / (aimf + af)
1608 * (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.))
1609 * (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.));
1611 G4double JPRFHEAVY = jprf0 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(af, 5.0 / 3.0) / IINERTTOT;
1612 G4double JPRFLIGHT = jprf0 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(aimf, 5.0 / 3.0) / IINERTTOT;
1613 if (af < 2.0) std::cout <<
"RN117-4,AF,ZF,EE,JPRFheavy" << std::endl;
1615 G4double vx1ev_imf = 0., vy1ev_imf = 0., vz1ev_imf = 0., zdummy = 0., adummy = 0.,
1616 tkedummy = 0., jprf1 = 0.;
1618 evapora(zf, af, &EEIMFP, JPRFHEAVY, &zff, &aff, &mtota, &vz1ev_imf, &vx1ev_imf, &vy1ev_imf,
1619 &FF11, &FIMF11, &zdummy, &adummy, &tkedummy, &jprf1, &inttype, &inum, EV_TEMP,
1620 &IEV_TAB_TEMP, &NbLamH);
1622 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1623 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1624 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1625 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1631 lorentz_boost(V_CM[0], V_CM[1], V_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4],
1632 &VXOUT, &VYOUT, &VZOUT);
1633 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1635 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
1636 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
1637 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
1639 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1644 opt->optimfallowed = 0;
1648 G4double zffimf, affimf, zdummy1 = 0., adummy1 = 0., tkedummy1 = 0., jprf2, vx2ev_imf,
1649 vy2ev_imf, vz2ev_imf;
1651 evapora(zimf, aimf, &EEIMF, JPRFLIGHT, &zffimf, &affimf, &mtota, &vz2ev_imf, &vx2ev_imf,
1652 &vy2ev_imf, &FF22, &FIMF22, &zdummy1, &adummy1, &tkedummy1, &jprf2, &inttype, &inum,
1653 EV_TEMP, &IEV_TAB_TEMP, &NbLamimf);
1655 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1656 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1657 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1658 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1664 lorentz_boost(V_CM[0], V_CM[1], V_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4],
1665 &VXOUT, &VYOUT, &VZOUT);
1666 lorentz_boost(VX2_IMF, VY2_IMF, VZ2_IMF, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1667 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
1668 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
1669 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
1671 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1684 lorentz_boost(VX2_IMF, VY2_IMF, VZ2_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1685 lorentz_boost(vx2ev_imf, vy2ev_imf, vz2ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1686 VIMF_CM[0] = VX2OUT;
1687 VIMF_CM[1] = VY2OUT;
1688 VIMF_CM[2] = VZ2OUT;
1693 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1694 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1695 VFP1_CM[0] = VX2OUT;
1696 VFP1_CM[1] = VY2OUT;
1697 VFP1_CM[2] = VZ2OUT;
1699 if (FF11 == 0 && FIMF11 == 0) {
1711 for (
G4int I = 0; I < 3; I++)
1714 else if (FF11 == 1 && FIMF11 == 0) {
1717 if (NbLam0 > 0) varntp->kfis = 20;
1719 opt->optimfallowed = 0;
1728 G4int IEV_TAB_FIS = 0, imode = 0;
1730 G4double vx1_fission = 0., vy1_fission = 0., vz1_fission = 0.;
1731 G4double vx2_fission = 0., vy2_fission = 0., vz2_fission = 0.;
1732 G4double vx_eva_sc = 0., vy_eva_sc = 0., vz_eva_sc = 0.;
1734 fission(af, zf, ee, jprf1, &vx1_fission, &vy1_fission, &vz1_fission, &vx2_fission,
1735 &vy2_fission, &vz2_fission, &ZFP1, &AFP1, &SFP1, &ZFP2, &AFP2, &SFP2, &imode,
1736 &vx_eva_sc, &vy_eva_sc, &vz_eva_sc, EV_TEMP, &IEV_TAB_FIS, &NbLamH);
1738 for (
G4int IJ = 0; IJ < IEV_TAB_FIS; IJ++) {
1739 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1740 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1741 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1747 lorentz_boost(VFP1_CM[0], VFP1_CM[1], VFP1_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3],
1748 EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1749 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1750 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1751 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1753 IEV_TAB = IEV_TAB + IEV_TAB_FIS;
1765 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1766 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1768 lorentz_boost(vx1_fission, vy1_fission, vz1_fission, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT,
1770 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1772 VFP1_CM[0] = VX2OUT;
1773 VFP1_CM[1] = VY2OUT;
1774 VFP1_CM[2] = VZ2OUT;
1783 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1784 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1786 lorentz_boost(vx2_fission, vy2_fission, vz2_fission, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT,
1788 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1790 VFP2_CM[0] = VX2OUT;
1791 VFP2_CM[1] = VY2OUT;
1792 VFP2_CM[2] = VZ2OUT;
1794 else if (FF11 == 0 && FIMF11 == 1) {
1797 opt->optimfallowed = 0;
1811 G4int NbLamimf1 = 0;
1812 G4double pbH1 = (af - zf) / (af - zf + aimf - zimf);
1813 for (
G4int i = 0; i < NbLamH; i++) {
1823 EkinR1 = tkeimf * aimf / (af + aimf);
1824 EkinR2 = tkeimf * af / (af + aimf);
1825 V1 = std::sqrt(EkinR1 / af) * 1.3887;
1826 V2 = std::sqrt(EkinR2 / aimf) * 1.3887;
1828 VPERP1 = std::sqrt(
V1 *
V1 - VZ1_IMFS * VZ1_IMFS);
1830 G4double VX1_IMFS = VPERP1 * std::sin(ALPHA1);
1831 G4double VY1_IMFS = VPERP1 * std::cos(ALPHA1);
1836 EEIMFP = ee * af / (af + aimf);
1837 EEIMF = ee * aimf / (af + aimf);
1841 0.40 * 931.490 * 1.160 * 1.160 * (std::pow(aimf, 5.0 / 3.0) + std::pow(af, 5.0 / 3.0))
1842 + 931.490 * 1.160 * 1.160 * aimf * af / (aimf + af)
1843 * (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.))
1844 * (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.));
1846 JPRFHEAVY = jprf1 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(af, 5.0 / 3.0) / IINERTTOT;
1847 JPRFLIGHT = jprf1 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(aimf, 5.0 / 3.0) / IINERTTOT;
1849 G4double zffs = 0., affs = 0., vx1ev_imfs = 0., vy1ev_imfs = 0., vz1ev_imfs = 0., jprf3 = 0.;
1851 evapora(zf, af, &EEIMFP, JPRFHEAVY, &zffs, &affs, &mtota, &vz1ev_imfs, &vx1ev_imfs,
1852 &vy1ev_imfs, &FF11, &FIMF11, &zdummy, &adummy, &tkedummy, &jprf3, &inttype, &inum,
1853 EV_TEMP, &IEV_TAB_TEMP, &NbLamH1);
1855 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1856 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1857 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1858 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1864 lorentz_boost(VFP1_CM[0], VFP1_CM[1], VFP1_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3],
1865 EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1866 lorentz_boost(vx1ev_imfs, vy1ev_imfs, vz1ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1868 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
1869 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
1870 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
1872 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1875 opt->optimfallowed = 0;
1881 G4double zffimfs = 0., affimfs = 0., vx2ev_imfs = 0., vy2ev_imfs = 0., vz2ev_imfs = 0.,
1884 evapora(zimf, aimf, &EEIMF, JPRFLIGHT, &zffimfs, &affimfs, &mtota, &vz2ev_imfs, &vx2ev_imfs,
1885 &vy2ev_imfs, &FF22, &FIMF22, &zdummy1, &adummy1, &tkedummy1, &jprf4, &inttype, &inum,
1886 EV_TEMP, &IEV_TAB_TEMP, &NbLamimf1);
1888 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1889 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1890 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1891 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1897 lorentz_boost(VFP1_CM[0], VFP1_CM[1], VFP1_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3],
1898 EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1899 lorentz_boost(vx2ev_imfs, vy2ev_imfs, vz2ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1901 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
1902 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
1903 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
1905 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1919 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1920 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1922 lorentz_boost(VX1_IMFS, VY1_IMFS, VZ1_IMFS, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT, &VZOUT);
1923 lorentz_boost(vx1ev_imfs, vy1ev_imfs, vz1ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1925 VFP1_CM[0] = VX2OUT;
1926 VFP1_CM[1] = VY2OUT;
1927 VFP1_CM[2] = VZ2OUT;
1934 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1935 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1937 lorentz_boost(VX2_IMFS, VY2_IMFS, VZ2_IMFS, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT, &VZOUT);
1938 lorentz_boost(vx2ev_imfs, vy2ev_imfs, vz2ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1940 VFP2_CM[0] = VX2OUT;
1941 VFP2_CM[1] = VY2OUT;
1942 VFP2_CM[2] = VZ2OUT;
1947 if (ftype != 1 && ftype != 21) {
1951 unstable_nuclei(AFP1, ZFP1, &afpnew, &zfpnew, IOUNSTABLE, VFP1_CM[0], VFP1_CM[1], VFP1_CM[2],
1952 &VP1X, &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
1954 if (IOUNSTABLE == 1) {
1960 for (
G4int I = 0; I < ILOOP; I++) {
1961 for (
G4int IJ = 0; IJ < 5; IJ++)
1962 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1964 IEV_TAB = IEV_TAB + ILOOP;
1970 unstable_nuclei(AFPIMF, ZFPIMF, &afpnew, &zfpnew, IOUNSTABLE, VIMF_CM[0], VIMF_CM[1],
1971 VIMF_CM[2], &VP1X, &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
1973 if (IOUNSTABLE == 1) {
1979 for (
G4int I = 0; I < ILOOP; I++) {
1980 for (
G4int IJ = 0; IJ < 5; IJ++)
1981 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1983 IEV_TAB = IEV_TAB + ILOOP;
1989 unstable_nuclei(AFP2, ZFP2, &afpnew, &zfpnew, IOUNSTABLE, VFP2_CM[0], VFP2_CM[1],
1990 VFP2_CM[2], &VP1X, &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
1992 if (IOUNSTABLE == 1) {
1998 for (
G4int I = 0; I < ILOOP; I++) {
1999 for (
G4int IJ = 0; IJ < 5; IJ++)
2000 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2002 IEV_TAB = IEV_TAB + ILOOP;
2009 if (ftype == 1 || ftype == 21) {
2013 unstable_nuclei(AFP1, ZFP1, &afpnew, &zfpnew, IOUNSTABLE, VFP1_CM[0], VFP1_CM[1], VFP1_CM[2],
2014 &VP1X, &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
2016 if (IOUNSTABLE == 1) {
2022 for (
G4int I = 0; I < ILOOP; I++) {
2023 for (
G4int IJ = 0; IJ < 5; IJ++)
2024 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2026 IEV_TAB = IEV_TAB + ILOOP;
2031 unstable_nuclei(AFP2, ZFP2, &afpnew, &zfpnew, IOUNSTABLE, VFP2_CM[0], VFP2_CM[1], VFP2_CM[2],
2032 &VP1X, &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
2034 if (IOUNSTABLE == 1) {
2040 for (
G4int I = 0; I < ILOOP; I++) {
2041 for (
G4int IJ = 0; IJ < 5; IJ++)
2042 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2044 IEV_TAB = IEV_TAB + ILOOP;
2050 unstable_nuclei(AFPIMF, ZFPIMF, &afpnew, &zfpnew, IOUNSTABLE, VIMF_CM[0], VIMF_CM[1],
2051 VIMF_CM[2], &VP1X, &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
2053 if (IOUNSTABLE == 1) {
2059 for (
G4int I = 0; I < ILOOP; I++) {
2060 for (
G4int IJ = 0; IJ < 5; IJ++)
2061 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2063 IEV_TAB = IEV_TAB + ILOOP;
2069 if ((ftype == 1 || ftype == 21) && (AFP2 <= 0 || AFP1 <= 0 || ZFP2 <= 0 || ZFP1 <= 0)) {
2070 std::cout <<
"ZFP1:" << ZFP1 << std::endl;
2071 std::cout <<
"AFP1:" << AFP1 << std::endl;
2072 std::cout <<
"ZFP2:" << ZFP2 << std::endl;
2073 std::cout <<
"AFP2:" << AFP2 << std::endl;
2077 EV_TAB[IEV_TAB][0] = ZFP1;
2078 EV_TAB[IEV_TAB][1] = AFP1;
2079 EV_TAB[IEV_TAB][5] = SFP1;
2080 EV_TAB[IEV_TAB][2] = VFP1_CM[0];
2081 EV_TAB[IEV_TAB][3] = VFP1_CM[1];
2082 EV_TAB[IEV_TAB][4] = VFP1_CM[2];
2083 IEV_TAB = IEV_TAB + 1;
2086 EV_TAB[IEV_TAB][0] = ZFP2;
2087 EV_TAB[IEV_TAB][1] = AFP2;
2088 EV_TAB[IEV_TAB][5] = SFP2;
2089 EV_TAB[IEV_TAB][2] = VFP2_CM[0];
2090 EV_TAB[IEV_TAB][3] = VFP2_CM[1];
2091 EV_TAB[IEV_TAB][4] = VFP2_CM[2];
2092 IEV_TAB = IEV_TAB + 1;
2096 EV_TAB[IEV_TAB][0] = ZFPIMF;
2097 EV_TAB[IEV_TAB][1] = AFPIMF;
2098 EV_TAB[IEV_TAB][5] = SFPIMF;
2099 EV_TAB[IEV_TAB][2] = VIMF_CM[0];
2100 EV_TAB[IEV_TAB][3] = VIMF_CM[1];
2101 EV_TAB[IEV_TAB][4] = VIMF_CM[2];
2102 IEV_TAB = IEV_TAB + 1;
3126 G4double problamb0 = (*problamb0_par);
3278 G4int choice_fisspart = 0;
3356 G4int fimf_allowed = opt->optimfallowed;
3375 mglms(a, zprf, fiss->optshp, &maz);
3376 mglms(a - 1.0, zprf, fiss->optshp, &ma1z);
3377 mglms(a - 2.0, zprf, fiss->optshp, &ma2z);
3378 mglms(a - 1.0, zprf - 1.0, fiss->optshp, &ma1z1);
3379 mglms(a - 2.0, zprf - 1.0, fiss->optshp, &ma2z1);
3380 mglms(a - 3.0, zprf - 1.0, fiss->optshp, &ma3z1);
3381 mglms(a - 3.0, zprf - 2.0, fiss->optshp, &ma3z2);
3382 mglms(a - 4.0, zprf - 2.0, fiss->optshp, &ma4z2);
3385 mglw(a, zprf, &maz);
3386 mglw(a - 1.0, zprf, &ma1z);
3387 mglw(a - 1.0, zprf - 1.0, &ma1z1);
3388 mglw(a - 2.0, zprf - 1.0, &ma2z1);
3389 mglw(a - 3.0, zprf - 1.0, &ma3z1);
3390 mglw(a - 3.0, zprf - 2.0, &ma3z2);
3391 mglw(a - 4.0, zprf - 2.0, &ma4z2);
3394 if ((a - 1.) == 3.0 && (zprf - 1.0) == 2.0) ma1z1 = -7.7181660;
3395 if ((a - 1.) == 4.0 && (zprf - 1.0) == 2.0) ma1z1 = -28.295992;
3400 sd = ma2z1 - maz - 2.2246;
3401 st = ma3z1 - maz - 8.481977;
3402 she = ma3z2 - maz - 7.7181660;
3403 sa = ma4z2 - maz - 28.295992;
3433 if (zprf <= 1.0e0 || a <= 1.0e0 || (a - zprf) < 0.0) {
3444 if (zprf <= 1.0e0 || a <= 2.0e0 || (a - zprf) < 1.0) {
3455 if (zprf <= 1.0e0 || a <= 3.0e0 || (a - zprf) < 2.0) {
3466 if (a - 4.0 <= 0.0 || zprf <= 2.0 || (a - zprf) < 2.0) {
3477 if (a - 3.0 <= 0.0 || zprf <= 2.0 || (a - zprf) < 1.0) {
3483 bhe =
max(bhe, 0.1);
3492 unbound(sn, sp, sd, st, she, sa, bp, bd, bt, bhe, ba, &probf, &probn, &probp, &probd, &probt,
3493 &probhe, &proba, &probimf, &probg, &ecn, &ecp, &ecd, &ect, &eche, &eca);
3499 if (fiss->ifis > 0) {
3504 barfit(k, k + j, il, &sbfis, &segs, &selmax);
3510 if (ef < 0.0) ef = 0.0;
3516 ef = ef + 0.51 * (1115. - 938. + sn - slamb0) / std::pow(a, 2. / 3.);
3525 xx =
fissility((k + j), k, NbLam0, sn, slamb0, fiss->optxfis);
3527 if (y < 0.0) y = 0.0;
3528 if (y > 1.0) y = 1.0;
3546 bshell = ecld->ecgnz[in][iz] - ecld->vgsld[in][iz];
3547 defbet = ecld->beta2[in][iz];
3549 iinert = 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a, 5.0 / 3.0)
3550 * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
3551 erot = jprf * jprf * 197.328 * 197.328 / (2. * iinert);
3554 bsbkbc(a, zprf, &bscn, &bkcn, &bccn);
3557 densniv(a, zprf, ee, 0.0, &densg, bshell, bscn, bkcn, &temp, fiss->optshp, fiss->optcol, defbet,
3558 &ecor, jprf, 0, &qrcn);
3610 lorb(a, a - 1., jprf, ee - sn, &dlout, &sdlout);
3611 djprf =
gausshaz(1, dlout, sdlout);
3612 if (IDjprf == 1) djprf = 0.0;
3613 jprfn = jprf + djprf;
3614 jprfn =
dint(std::abs(jprfn));
3616 bshell = ecld->ecgnz[ind][izd] - ecld->vgsld[ind][izd];
3617 defbet = ecld->beta2[ind][izd];
3619 iinert = 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a - 1., 5.0 / 3.0)
3620 * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
3621 erotn = jprfn * jprfn * 197.328 * 197.328 / (2. * iinert);
3622 bsbkbc(a - 1., zprf, &bs, &bk, &bc);
3625 densniv(a - 1.0, zprf, ee, sn, &densn, bshell, bs, bk, &temp, fiss->optshp, fiss->optcol,
3626 defbet, &ecor, jprfn, 0, &qr);
3631 if (imaxwell == 1) {
3637 std::cout <<
"WARNING: FVMAXHAZ_NEUT CALLED MORE THAN 100 TIMES" << std::endl;
3640 if (ecn > (ee - sn)) {
3641 if ((ee - sn) < rnt)
3646 if (ecn <= 0.0)
goto dir1234;
3668 lorb(a, a - 1., jprf, ee - sbp, &dlout, &sdlout);
3669 djprf =
gausshaz(1, dlout, sdlout);
3670 if (IDjprf == 1) djprf = 0.0;
3671 jprfp = jprf + djprf;
3672 jprfp =
dint(std::abs(jprfp));
3674 bshell = ecld->ecgnz[ind][izd] - ecld->vgsld[ind][izd];
3675 defbet = ecld->beta2[ind][izd];
3677 iinert = 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a - 1., 5.0 / 3.0)
3678 * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
3679 erotp = jprfp * jprfp * 197.328 * 197.328 / (2. * iinert);
3681 bsbkbc(a - 1., zprf - 1., &bs, &bk, &bc);
3684 densniv(a - 1.0, zprf - 1.0, ee, sbp, &densp, bshell, bs, bk, &temp, fiss->optshp, fiss->optcol,
3685 defbet, &ecor, jprfp, 0, &qr);
3690 if (imaxwell == 1) {
3696 std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3699 if (ecp > (ee - sbp)) {
3700 if ((ee - sbp) < rpt)
3705 if (ecp <= 0.0)
goto dir1235;
3709 ecp = 2.0 * pt + bp;
3724 if ((in >= 2) && (iz >= 2)) {
3728 lorb(a, a - 2., jprf, ee - sbd, &dlout, &sdlout);
3729 djprf =
gausshaz(1, dlout, sdlout);
3730 if (IDjprf == 1) djprf = 0.0;
3731 jprfd = jprf + djprf;
3732 jprfd =
dint(std::abs(jprfd));
3734 bshell = ecld->ecgnz[ind][izd] - ecld->vgsld[ind][izd];
3735 defbet = ecld->beta2[ind][izd];
3737 iinert = 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a - 2., 5.0 / 3.0)
3738 * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
3739 erotd = jprfd * jprfd * 197.328 * 197.328 / (2. * iinert);
3741 bsbkbc(a - 2., zprf - 1., &bs, &bk, &bc);
3744 densniv(a - 2.0, zprf - 1.0e0, ee, sbd, &densd, bshell, bs, bk, &temp, fiss->optshp,
3745 fiss->optcol, defbet, &ecor, jprfd, 0, &qr);
3751 if (imaxwell == 1) {
3757 std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3760 if (ecd > (ee - sbd)) {
3761 if ((ee - sbd) < rdt)
3766 if (ecd <= 0.0)
goto dir1236;
3770 ecd = 2.0 * dt + bd;
3785 if ((in >= 3) && (iz >= 2)) {
3789 lorb(a, a - 3., jprf, ee - sbt, &dlout, &sdlout);
3790 djprf =
gausshaz(1, dlout, sdlout);
3791 if (IDjprf == 1) djprf = 0.0;
3792 jprft = jprf + djprf;
3793 jprft =
dint(std::abs(jprft));
3795 bshell = ecld->ecgnz[ind][izd] - ecld->vgsld[ind][izd];
3796 defbet = ecld->beta2[ind][izd];
3798 iinert = 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a - 3., 5.0 / 3.0)
3799 * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
3800 erott = jprft * jprft * 197.328 * 197.328 / (2. * iinert);
3802 bsbkbc(a - 3., zprf - 1., &bs, &bk, &bc);
3805 densniv(a - 3.0, zprf - 1.0, ee, sbt, &denst, bshell, bs, bk, &temp, fiss->optshp, fiss->optcol,
3806 defbet, &ecor, jprft, 0, &qr);
3812 if (imaxwell == 1) {
3818 std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3821 if (ect > (ee - sbt)) {
3822 if ((ee - sbt) < rtt)
3827 if (ect <= 0.0)
goto dir1237;
3831 ect = 2.0 * tt + bt;
3846 if ((in >= 3) && (iz >= 3)) {
3850 lorb(a, a - 4., jprf, ee - sba, &dlout, &sdlout);
3851 djprf =
gausshaz(1, dlout, sdlout);
3852 if (IDjprf == 1) djprf = 0.0;
3853 jprfa = jprf + djprf;
3854 jprfa =
dint(std::abs(jprfa));
3856 bshell = ecld->ecgnz[ind][izd] - ecld->vgsld[ind][izd];
3857 defbet = ecld->beta2[ind][izd];
3859 iinert = 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a - 4., 5.0 / 3.0)
3860 * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
3861 erota = jprfa * jprfa * 197.328 * 197.328 / (2. * iinert);
3863 bsbkbc(a - 4., zprf - 2., &bs, &bk, &bc);
3866 densniv(a - 4.0, zprf - 2.0, ee, sba, &densa, bshell, bs, bk, &temp, fiss->optshp, fiss->optcol,
3867 defbet, &ecor, jprfa, 0, &qr);
3873 if (imaxwell == 1) {
3879 std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3882 if (eca > (ee - sba)) {
3883 if ((ee - sba) < rat)
3888 if (eca <= 0.0)
goto dir1238;
3892 eca = 2.0 * at + ba;
3907 if ((in >= 2) && (iz >= 3)) {
3911 lorb(a, a - 3., jprf, ee - sbhe, &dlout, &sdlout);
3912 djprf =
gausshaz(1, dlout, sdlout);
3913 if (IDjprf == 1) djprf = 0.0;
3914 jprfhe = jprf + djprf;
3915 jprfhe =
dint(std::abs(jprfhe));
3917 bshell = ecld->ecgnz[ind][izd] - ecld->vgsld[ind][izd];
3918 defbet = ecld->beta2[ind][izd];
3920 iinert = 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a - 3., 5.0 / 3.0)
3921 * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
3922 erothe = jprfhe * jprfhe * 197.328 * 197.328 / (2. * iinert);
3924 bsbkbc(a - 3., zprf - 2., &bs, &bk, &bc);
3927 densniv(a - 3.0, zprf - 2.0, ee, sbhe, &denshe, bshell, bs, bk, &temp, fiss->optshp,
3928 fiss->optcol, defbet, &ecor, jprfhe, 0, &qr);
3934 if (imaxwell == 1) {
3940 std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3943 if (eche > (ee - sbhe)) {
3944 if ((ee - sbhe) < rhet)
3949 if (eche <= 0.0)
goto dir1239;
3953 eche = 2.0 * het + bhe;
3969 if (in >= 2 && NbLam0 > 0) {
3973 lorb(a, a - 1., jprf, ee - slamb0, &dlout, &sdlout);
3974 djprf =
gausshaz(1, dlout, sdlout);
3975 if (IDjprf == 1) djprf = 0.0;
3976 jprflamb0 = jprf + djprf;
3977 jprflamb0 =
dint(std::abs(jprflamb0));
3979 bshell = ecld->ecgnz[ind][izd] - ecld->vgsld[ind][izd];
3980 defbet = ecld->beta2[ind][izd];
3982 iinert = 0.4 * 931.49 * 1.16 * 1.16 * std::pow(a - 1., 5.0 / 3.0)
3983 * (1.0 + 0.5 * std::sqrt(5. / (4. * pi)) * defbet);
3984 erotlamb0 = jprflamb0 * jprflamb0 * 197.328 * 197.328 / (2. * iinert);
3985 bsbkbc(a - 1., zprf, &bs, &bk, &bc);
3988 densniv(a - 1.0, zprf, ee, slamb0, &denslamb0, bshell, bs, bk, &temp, fiss->optshp,
3989 fiss->optcol, defbet, &ecor, jprflamb0, 0, &qr);
3992 if (denslamb0 > 0.) {
3994 if (imaxwell == 1) {
4000 std::cout <<
"WARNING: FVMAXHAZ_NEUT CALLED MORE THAN 100 TIMES" << std::endl;
4003 if (eclamb0 > (ee - slamb0)) {
4004 if ((ee - slamb0) < rlamb0t)
4005 eclamb0 = ee - slamb0;
4009 if (eclamb0 <= 0.0)
goto dir1240;
4012 eclamb0 = 2.0 * lamb0t;
4040 gn =
width(a, zprf, 1.0, 0.0, nt, 0.0, sn, ee - erotn) * densn / densg;
4047 width(a, zprf, 1.0, 1.0, pt, bp, sbp, ee - erotp) * densp / densg *
pen(a, 1.0, omegap, pt);
4054 width(a, zprf, 2.0, 1.0, dt, bd, sbd, ee - erotd) * densd / densg *
pen(a, 2.0, omegad, dt);
4061 width(a, zprf, 3.0, 1.0, tt, bt, sbt, ee - erott) * denst / densg *
pen(a, 3.0, omegat, tt);
4063 if (denshe <= 0.0) {
4067 ghe =
width(a, zprf, 3.0, 2.0, het, bhe, sbhe, ee - erothe) * denshe / densg
4068 *
pen(a, 3.0, omegahe, het);
4075 width(a, zprf, 4.0, 2.0, at, ba, sba, ee - erota) * densa / densg *
pen(a, 4.0, omegaa, at);
4077 if (denslamb0 <= 0.0) {
4081 glamb0 =
width(a, zprf, 1.0, -2.0, lamb0t, 0.0, slamb0, ee - erotlamb0) * denslamb0 / densg;
4089 G4int izcn = 0, incn = 0, inmin = 0, inmax = 0, inmi = 0, inma = 0;
4092 if (fimf_allowed == 0 || zprf <= 5.0 || a <= 7.0) {
4099 mglms(a, zprf, opt->optshpimf, &mazz);
4113 inmin =
max(inmin, incn - inma);
4114 inmax =
min(inmax, incn - inmi);
4116 inmax =
max(inmax, inmin);
4118 for (
G4int iaimf = izimf + inmin; iaimf <= izimf + inmax; iaimf++) {
4120 if (aimf >= a || zimf >= zprf) {
4125 mglms(a - aimf, zprf - zimf, opt->optshpimf, &mares);
4126 mglms(aimf, zimf, opt->optshpimf, &maimf);
4131 defbetimf = ecld->beta2[
idnint(aimf - zimf)][
idnint(zimf)]
4132 + ecld->beta2[
idnint(a - aimf - zprf + zimf)][
idnint(zprf - zimf)];
4134 iinert = 0.40 * 931.490 * 1.160 * 1.160 * std::pow(a, 5.0 / 3.0)
4135 * (std::pow(aimf, 5.0 / 3.0) + std::pow(a - aimf, 5.0 / 3.0))
4136 + 931.490 * 1.160 * 1.160 * aimf * (a - aimf) / a
4137 * (std::pow(aimf, 1.0 / 3.0) + std::pow(a - aimf, 1.0 / 3.0))
4138 * (std::pow(aimf, 1.0 / 3.0) + std::pow(a - aimf, 1.0 / 3.0));
4140 erot = jprf * jprf * 197.328 * 197.328 / (2.0 * iinert);
4143 if (densg == 0.0 || ee < (sbimf + erot)) {
4151 densniv(a, zprf, ee, sbimf, &densimf, 0.0, bsimf, 1.0, &timf, 0, 0, defbetimf, &ecor,
4154 imfarg = (sbimf + erotcn - erot) / timf;
4155 if (imfarg > 200.0) imfarg = 200.0;
4166 width_imf =
width(a, zprf, aimf, zimf, timf, bimf, sbimf, ee - erot) * std::exp(-imfarg)
4170 gimf3 = gimf3 + width_imf;
4186 inmin =
max(inmin, incn - inma);
4187 inmax =
min(inmax, incn - inmi);
4189 inmax =
max(inmax, inmin);
4191 for (
G4int iaimf = izimf + inmin; iaimf <= izimf + inmax; iaimf++) {
4193 if (aimf >= a || zimf >= zprf) {
4198 mglms(a - aimf, zprf - zimf, opt->optshpimf, &mares);
4199 mglms(aimf, zimf, opt->optshpimf, &maimf);
4204 defbetimf = ecld->beta2[
idnint(aimf - zimf)][
idnint(zimf)]
4205 + ecld->beta2[
idnint(a - aimf - zprf + zimf)][
idnint(zprf - zimf)];
4207 iinert = 0.40 * 931.490 * 1.160 * 1.160 * std::pow(a, 5.0 / 3.0)
4208 * (std::pow(aimf, 5.0 / 3.0) + std::pow(a - aimf, 5.0 / 3.0))
4209 + 931.490 * 1.160 * 1.160 * aimf * (a - aimf) / a
4210 * (std::pow(aimf, 1.0 / 3.0) + std::pow(a - aimf, 1.0 / 3.0))
4211 * (std::pow(aimf, 1.0 / 3.0) + std::pow(a - aimf, 1.0 / 3.0));
4213 erot = jprf * jprf * 197.328 * 197.328 / (2.0 * iinert);
4216 if (densg == 0.0 || ee < (sbimf + erot)) {
4224 densniv(a, zprf, ee, sbimf, &densimf, 0.0, bsimf, 1.0, &timf, 0, 0, defbetimf, &ecor,
4227 imfarg = (sbimf + erotcn - erot) / timf;
4228 if (imfarg > 200.0) imfarg = 200.0;
4238 width_imf =
width(a, zprf, aimf, zimf, timf, bimf, sbimf, ee - erot) * std::exp(-imfarg)
4242 gimf5 = gimf5 + width_imf;
4247 if (gimf3 <= 0.0 || gimf5 <= 0.0) {
4254 b_imf = (std::log10(gimf3) - std::log10(gimf5)) / (std::log10(3.0) - std::log10(5.0));
4256 if (b_imf >= -1.01) b_imf = -1.01;
4257 if (b_imf <= -100.0) {
4264 a_imf = gimf3 / std::pow(3.0, b_imf);
4265 gimf = a_imf * (std::pow(zprf, b_imf + 1.0) - std::pow(3.0, b_imf + 1.0)) / (b_imf + 1.0);
4269 if (gimf < 1.e-10) gimf = 0.0;
4276 pa = (ald->av) * a + (ald->as) * std::pow(a, 2. / 3.) + (ald->ak) * std::pow(a, 1. / 3.);
4277 gamma = 2.5 * pa * std::pow(a, -4. / 3.);
4278 gfactor = 1. + gamma * ecld->ecgnz[in][iz];
4279 if (gfactor <= 0.) {
4283 gtemp = 17.60 / (std::pow(a, 0.699) * std::sqrt(gfactor));
4287 gg = 0.624e-9 * std::pow(a, 1.6) * std::pow(gtemp, 5.);
4291 if (gammaemission == 1) {
4297 gsum = ga + ghe + gd + gt + gp + gn + gimf + gg + glamb0;
4310 if (fiss->ifis == 0 || (zprf * zprf / a <= 22.74 && zprf < 60.)) {
4318 fission_width(zprf, a, ee, bssp, bksp, ef, y, &gf, &temp, jprf, 0, 1, fiss->optcol,
4319 fiss->optshp, densg);
4338 if (fiss->bet <= 0.) {
4339 gtotal = ga + ghe + gp + gd + gt + gn + gg + gimf + gf + glamb0;
4340 if (gtotal <= 0.0) {
4354 probf = gf / gtotal;
4355 probn = gn / gtotal;
4356 probp = gp / gtotal;
4357 probd = gd / gtotal;
4358 probt = gt / gtotal;
4359 probhe = ghe / gtotal;
4360 proba = ga / gtotal;
4361 probg = gg / gtotal;
4362 probimf = gimf / gtotal;
4363 problamb0 = glamb0 / gtotal;
4377 fomega_sp(a, y, &mfcd, &omegasp, &homegasp);
4378 cf =
cram((NbLam0 > 0 ? fiss->bethyp : fiss->bet), homegasp);
4381 fomega_gs(a, zprf, &k1, &omegags, &homegags);
4382 tauc =
tau((NbLam0 > 0 ? fiss->bethyp : fiss->bet), homegags, ef, ft);
4395 part_fiss((NbLam0 > 0 ? fiss->bethyp : fiss->bet), gsum, gf, y, tauc, ts1, tsum, &choice_fisspart,
4396 zprf, a, ft, &t_lapse, &gff);
4401 tsum = tsum + t_lapse;
4404 if (choice_fisspart == 2) {
4422 gtotal = ga + ghe + gp + gd + gt + gn + gimf + gg + glamb0;
4423 if (gtotal <= 0.0) {
4438 probn = gn / gtotal;
4439 probp = gp / gtotal;
4440 probd = gd / gtotal;
4441 probt = gt / gtotal;
4442 probhe = ghe / gtotal;
4443 proba = ga / gtotal;
4444 probg = gg / gtotal;
4445 probimf = gimf / gtotal;
4446 problamb0 = glamb0 / gtotal;
4452 gtotal = ga + ghe + gp + gd + gt + gn + gg + gimf + glamb0;
4453 if (gtotal <= 0.0) {
4467 probn = gn / gtotal;
4468 probp = gp / gtotal;
4469 probd = gd / gtotal;
4470 probt = gt / gtotal;
4471 probhe = ghe / gtotal;
4472 proba = ga / gtotal;
4473 probg = gg / gtotal;
4474 probimf = gimf / gtotal;
4475 problamb0 = glamb0 / gtotal;
4479 ptotl = probp + probd + probt + probn + probhe + proba + probg + probimf + probf + problamb0;
4485 (*probp_par) = probp;
4486 (*probd_par) = probd;
4487 (*probt_par) = probt;
4488 (*probn_par) = probn;
4489 (*probhe_par) = probhe;
4490 (*proba_par) = proba;
4491 (*probg_par) = probg;
4492 (*probimf_par) = probimf;
4493 (*problamb0_par) = problamb0;
4494 (*probf_par) = probf;
4495 (*ptotl_par) = ptotl;
4502 (*slamb0_par) = slamb0;
4515 (*eclamb0_par) = eclamb0;
4523 (*jprfn_par) = jprfn;
4524 (*jprfp_par) = jprfp;
4525 (*jprfd_par) = jprfd;
4526 (*jprft_par) = jprft;
4527 (*jprfhe_par) = jprfhe;
4528 (*jprfa_par) = jprfa;
4529 (*jprflamb0_par) = jprflamb0;
5591 for (
G4int init_i = 0; init_i < 7; init_i++) {
5595 for (
G4int init_i = 0; init_i < 10; init_i++) {
5599 G4double a = 0.0, z = 0.0, amin = 0.0, amax = 0.0, amin2 = 0.0;
5600 G4double amax2 = 0.0, aa = 0.0, zz = 0.0, bfis = 0.0;
5601 G4double bfis0 = 0.0, ell = 0.0, el = 0.0, egs = 0.0, el80 = 0.0, el20 = 0.0;
5602 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;
5603 G4double aj = 0.0, ak = 0.0, a1 = 0.0, a2 = 0.0;
5605 G4int i = 0, j = 0, k = 0, m = 0;
5610 G4double emncof[4][5] = {{-9.01100e+2, -1.40818e+3, 2.77000e+3, -7.06695e+2, 8.89867e+2},
5611 {1.35355e+4, -2.03847e+4, 1.09384e+4, -4.86297e+3, -6.18603e+2},
5612 {-3.26367e+3, 1.62447e+3, 1.36856e+3, 1.31731e+3, 1.53372e+2},
5613 {7.48863e+3, -1.21581e+4, 5.50281e+3, -1.33630e+3, 5.05367e-2}};
5615 G4double elmcof[4][5] = {{1.84542e+3, -5.64002e+3, 5.66730e+3, -3.15150e+3, 9.54160e+2},
5616 {-2.24577e+3, 8.56133e+3, -9.67348e+3, 5.81744e+3, -1.86997e+3},
5617 {2.79772e+3, -8.73073e+3, 9.19706e+3, -4.91900e+3, 1.37283e+3},
5618 {-3.01866e+1, 1.41161e+3, -2.85919e+3, 2.13016e+3, -6.49072e+2}};
5621 {9.43596e4, -2.241997e5, 2.223237e5, -1.324408e5, 4.68922e4, -8.83568e3},
5622 {-1.655827e5, 4.062365e5, -4.236128e5, 2.66837e5, -9.93242e4, 1.90644e4},
5623 {1.705447e5, -4.032e5, 3.970312e5, -2.313704e5, 7.81147e4, -1.322775e4},
5624 {-9.274555e4, 2.278093e5, -2.422225e5, 1.55431e5, -5.78742e4, 9.97505e3}};
5626 G4double elzcof[7][7] = {{5.11819909e+5, -1.30303186e+6, 1.90119870e+6, -1.20628242e+6,
5627 5.68208488e+5, 5.48346483e+4, -2.45883052e+4},
5628 {-1.13269453e+6, 2.97764590e+6, -4.54326326e+6, 3.00464870e+6,
5629 -1.44989274e+6, -1.02026610e+5, 6.27959815e+4},
5630 {1.37543304e+6, -3.65808988e+6, 5.47798999e+6, -3.78109283e+6,
5631 1.84131765e+6, 1.53669695e+4, -6.96817834e+4},
5632 {-8.56559835e+5, 2.48872266e+6, -4.07349128e+6, 3.12835899e+6,
5633 -1.62394090e+6, 1.19797378e+5, 4.25737058e+4},
5634 {3.28723311e+5, -1.09892175e+6, 2.03997269e+6, -1.77185718e+6,
5635 9.96051545e+5, -1.53305699e+5, -1.12982954e+4},
5636 {4.15850238e+4, 7.29653408e+4, -4.93776346e+5, 6.01254680e+5,
5637 -4.01308292e+5, 9.65968391e+4, -3.49596027e+3},
5638 {-1.82751044e+5, 3.91386300e+5, -3.03639248e+5, 1.15782417e+5,
5639 -4.24399280e+3, -6.11477247e+3, 3.66982647e+2}};
5641 const G4int sizex = 5;
5642 const G4int sizey = 6;
5643 const G4int sizez = 4;
5645 G4double egscof[sizey][sizey][sizez];
5647 G4double egs1[sizey][sizex] = {{1.927813e5, 7.666859e5, 6.628436e5, 1.586504e5, -7.786476e3},
5648 {-4.499687e5, -1.784644e6, -1.546968e6, -4.020658e5, -3.929522e3},
5649 {4.667741e5, 1.849838e6, 1.641313e6, 5.229787e5, 5.928137e4},
5650 {-3.017927e5, -1.206483e6, -1.124685e6, -4.478641e5, -8.682323e4},
5651 {1.226517e5, 5.015667e5, 5.032605e5, 2.404477e5, 5.603301e4},
5652 {-1.752824e4, -7.411621e4, -7.989019e4, -4.175486e4, -1.024194e4}};
5654 G4double egs2[sizey][sizex] = {{-6.459162e5, -2.903581e6, -3.048551e6, -1.004411e6, -6.558220e4},
5655 {1.469853e6, 6.564615e6, 6.843078e6, 2.280839e6, 1.802023e5},
5656 {-1.435116e6, -6.322470e6, -6.531834e6, -2.298744e6, -2.639612e5},
5657 {8.665296e5, 3.769159e6, 3.899685e6, 1.520520e6, 2.498728e5},
5658 {-3.302885e5, -1.429313e6, -1.512075e6, -6.744828e5, -1.398771e5},
5659 {4.958167e4, 2.178202e5, 2.400617e5, 1.167815e5, 2.663901e4}};
5661 G4double egs3[sizey][sizex] = {{3.117030e5, 1.195474e6, 9.036289e5, 6.876190e4, -6.814556e4},
5662 {-7.394913e5, -2.826468e6, -2.152757e6, -2.459553e5, 1.101414e5},
5663 {7.918994e5, 3.030439e6, 2.412611e6, 5.228065e5, 8.542465e3},
5664 {-5.421004e5, -2.102672e6, -1.813959e6, -6.251700e5, -1.184348e5},
5665 {2.370771e5, 9.459043e5, 9.026235e5, 4.116799e5, 1.001348e5},
5666 {-4.227664e4, -1.738756e5, -1.795906e5, -9.292141e4, -2.397528e4}};
5668 G4double egs4[sizey][sizex] = {{-1.072763e5, -5.973532e5, -6.151814e5, 7.371898e4, 1.255490e5},
5669 {2.298769e5, 1.265001e6, 1.252798e6, -2.306276e5, -2.845824e5},
5670 {-2.093664e5, -1.100874e6, -1.009313e6, 2.705945e5, 2.506562e5},
5671 {1.274613e5, 6.190307e5, 5.262822e5, -1.336039e5, -1.115865e5},
5672 {-5.715764e4, -2.560989e5, -2.228781e5, -3.222789e3, 1.575670e4},
5673 {1.189447e4, 5.161815e4, 4.870290e4, 1.266808e4, 2.069603e3}};
5675 for (i = 0; i < sizey; i++) {
5676 for (j = 0; j < sizex; j++) {
5677 egscof[i][j][0] = egs1[i][j];
5678 egscof[i][j][1] = egs2[i][j];
5679 egscof[i][j][2] = egs3[i][j];
5680 egscof[i][j][3] = egs4[i][j];
5685 if (iz < 19 || iz > 122) {
5689 if (iz > 122 && il > 0) {
5696 amin = 1.2e0 * z + 0.01e0 * z * z;
5697 amax = 5.8e0 * z - 0.024e0 * z * z;
5699 if (a < amin || a > amax) {
5711 for (i = 0; i < 7; i++) {
5712 for (j = 0; j < 7; j++) {
5713 bfis0 = bfis0 + elzcof[j][i] * pz[i] * pa[j];
5719 if ((fiss->optshp == 1) || (fiss->optshp == 3)) {
5720 bfis = bfis - ecld->ecgnz[in][iz];
5723 if (
mod(in, 2) == 1) {
5724 bfis = bfis * (4.5114 - 2.2687 * in / iz);
5727 bfis = bfis * (3.3931 - 1.5338 * in / iz);
5731 if (in * 1.0 / iz > 1.52) bfis = bfis * (1.1222 - 0.10886 * (in) / iz) - 0.1;
5733 if (iz >= 94 && iz <= 98 && in < 156) {
5735 if (
mod(in, 2) == 0 &&
mod(iz, 2) == 0) {
5737 bfis = bfis - (11.54108 * in / iz - 18.074);
5741 if (
mod(in, 2) == 1 &&
mod(iz, 2) == 1) {
5743 bfis = bfis - (14.567 * in / iz - 23.266);
5747 if (
mod(in, 2) == 0 &&
mod(iz, 2) == 1) {
5749 bfis = bfis - (13.662 * in / iz - 21.656);
5753 if (
mod(in, 2) == 1 &&
mod(iz, 2) == 0) {
5755 bfis = bfis - (13.662 * in / iz - 21.656);
5759 else if (iz == 84) {
5762 else if (iz == 83) {
5765 else if (iz == 82) {
5768 else if (iz == 81) {
5771 else if (iz == 80) {
5774 else if (iz <= 79) {
5785 amin2 = 1.4e0 * z + 0.009e0 * z * z;
5786 amax2 = 20.e0 + 3.0e0 * z;
5788 if ((a < amin2 - 5.e0 || a > amax2 + 10.e0) && il > 0) {
5798 for (i = 0; i < 4; i++) {
5799 for (j = 0; j < 5; j++) {
5800 el80 = el80 + elmcof[i][j] * pz[j] * pa[i];
5801 el20 = el20 + emncof[i][j] * pz[j] * pa[i];
5812 for (i = 0; i < 4; i++) {
5813 for (j = 0; j < 6; j++) {
5814 elmax = elmax + emxcof[i][j] * pz[j] * pa[i];
5825 x = sel20 / (*selmax);
5826 y = sel80 / (*selmax);
5830 q = 0.2 / (std::pow(sel20, 2) * std::pow(sel80, 2) * (sel20 - sel80));
5831 qa = q * (4.0 * std::pow(sel80, 3) - std::pow(sel20, 3));
5832 qb = -q * (4.0 * std::pow(sel80, 2) - std::pow(sel20, 2));
5833 bfis = bfis * (1.0 + qa * std::pow(el, 2) + qb * std::pow(el, 3));
5837 aj = (-20.0 * std::pow(x, 5) + 25.e0 * std::pow(x, 4) - 4.0) * std::pow((y - 1.0), 2) * y * y;
5838 ak = (-20.0 * std::pow(y, 5) + 25.0 * std::pow(y, 4) - 1.0) * std::pow((x - 1.0), 2) * x * x;
5839 q = 0.2 / (std::pow((y - x) * ((1.0 - x) * (1.0 - y) * x * y), 2));
5840 qa = q * (aj * y - ak * x);
5841 qb = -q * (aj * (2.0 * y + 1.0) - ak * (2.0 * x + 1.0));
5843 a1 = 4.0 * std::pow(z, 5) - 5.0 * std::pow(z, 4) + 1.0;
5844 a2 = qa * (2.e0 * z + 1.e0);
5845 bfis = bfis * (a1 + (z - 1.e0) * (a2 + qb * z) * z * z * (z - 1.e0));
5852 if (el > (*selmax)) {
5858 if (el > (*selmax)) {
5862 for (k = 0; k < 4; k++) {
5863 for (l = 0; l < 6; l++) {
5864 for (m = 0; m < 5; m++) {
5865 egs = egs + egscof[l][m][k] * pz[l] * pa[k] * pl[2 * m];
5871 if ((*segs) < 0.0) {
8290 Delta_U1_shell_max = -2.45;
8296 Delta_U2_shell = -2.45;
8306 G4double Fwidth_asymm1, Fwidth_asymm2, Fwidth_symm;
8308 Fwidth_asymm1 = 0.65;
8309 Fwidth_asymm2 = 0.65;
8351 Friction_factor = 1.0;
8354 G4double cN_asymm1_shell, cN_asymm2_shell;
8355 G4double gamma, gamma_heavy1, gamma_heavy2;
8371 G4double Sasymm1 = 0., Sasymm2 = 0., Ssymm = 0., Ysum = 0., Yasymm = 0.;
8373 G4double wNasymm1_saddle, wNasymm2_saddle, wNsymm_saddle;
8374 G4double wNasymm2_scission, wNsymm_scission;
8375 G4double wNasymm1, wNasymm2, wNsymm;
8389 G4double beta = 0., beta1 = 0., beta2 = 0.;
8398 G4double A_levdens_heavy1, A_levdens_heavy2;
8402 G4double epsilon_1_saddle, epsilon0_1_saddle;
8403 G4double epsilon_2_saddle, epsilon0_2_saddle, epsilon_symm_saddle;
8408 G4double E_eff1_saddle, E_eff2_saddle;
8409 G4double Epot0_mode1_saddle, Epot0_mode2_saddle, Epot0_symm_saddle;
8410 G4double Epot_mode1_saddle, Epot_mode2_saddle, Epot_symm_saddle;
8411 G4double E_defo, E_defo1, E_defo2, E_scission_pre = 0., E_scission_post;
8417 G4double MassCurv_scis, MassCurv_sadd;
8419 G4double Nheavy1_shell, Nheavy2_shell;
8424 G4double Z_scission, N_scission, A_scission;
8426 G4double beta1gs = 0., beta2gs = 0., betags = 0.;
8428 G4double DSN132, Delta_U1_shell, E_eff0_saddle;
8429 G4int NbLam0 = (*NbLam0_par);
8434 cN_asymm1_shell = 0.700 *
N / Z;
8435 cN_asymm2_shell = 0.040 *
N / Z;
8439 DSN132 = Nheavy1_in -
N / Z * Zheavy1_in;
8440 Aheavy1 = Nheavy1_in + Zheavy1_in + 0.340 * DSN132;
8449 Delta_U1_shell = Delta_U1_shell_max + U1NZ_SLOPE * std::abs(DSN132);
8450 Delta_U1_shell =
min(0., Delta_U1_shell);
8454 Nheavy1 =
N /
A * Aheavy1;
8455 Aheavy2 = Nheavy2 *
A /
N;
8459 A_levdens =
A / xLevdens;
8460 gamma = A_levdens / (0.40 * std::pow(
A, 1.3333)) * FGAMMA;
8461 A_levdens_heavy1 = Aheavy1 / xLevdens;
8462 gamma_heavy1 = A_levdens_heavy1 / (0.40 * std::pow(Aheavy1, 1.3333)) * FGAMMA * FGAMMA1;
8463 A_levdens_heavy2 = Aheavy2 / xLevdens;
8464 gamma_heavy2 = A_levdens_heavy2 / (0.40 * std::pow(Aheavy2, 1.3333)) * FGAMMA;
8468 E_saddle_scission = (-24. + 0.02227 * Z * Z / std::pow(
A, 0.33333)) * Friction_factor;
8469 E_saddle_scission =
max(0.0, E_saddle_scission);
8475 Z2_over_A_eff = Z * Z /
A;
8477 if (Z2_over_A_eff < 34.0)
8478 MassCurv_scis = std::pow(10., -1.093364 + 0.082933 * Z2_over_A_eff
8479 - 0.0002602 * Z2_over_A_eff * Z2_over_A_eff);
8481 MassCurv_scis = std::pow(10., 3.053536 - 0.056477 * Z2_over_A_eff
8482 + 0.0002454 * Z2_over_A_eff * Z2_over_A_eff);
8487 MassCurv_sadd = X_s2s * MassCurv_scis;
8489 cN_symm = 8.0 / std::pow(
N, 2.) * MassCurv_scis;
8490 cN_symm_sadd = 8.0 / std::pow(
N, 2.) * MassCurv_sadd;
8492 Nheavy1_shell = Nheavy1;
8496 (cN_symm_sadd * Nsymm
8497 + cN_asymm1_shell *
Uwash(E /
A * Aheavy1, Ecrit, FREDSHELL, gamma_heavy1) * Nheavy1_shell)
8498 / (cN_symm_sadd + cN_asymm1_shell *
Uwash(E /
A * Aheavy1, Ecrit, FREDSHELL, gamma_heavy1));
8501 (cN_symm_sadd * Nsymm + cN_asymm1_shell * Nheavy1_shell) / (cN_symm_sadd + cN_asymm1_shell);
8504 Nheavy2_NZ = Nheavy2;
8505 Nheavy2_shell = Nheavy2_NZ;
8508 (cN_symm_sadd * Nsymm
8509 + cN_asymm2_shell *
Uwash(E /
A * Aheavy2, Ecrit, FREDSHELL, gamma_heavy2) * Nheavy2_shell)
8510 / (cN_symm_sadd + cN_asymm2_shell *
Uwash(E /
A * Aheavy2, Ecrit, FREDSHELL, gamma_heavy2));
8513 (cN_symm_sadd * Nsymm + cN_asymm2_shell * Nheavy2_shell) / (cN_symm_sadd + cN_asymm2_shell);
8515 Delta_U1 = Delta_U1_shell
8516 + (Nheavy1_shell - Nheavy1_eff) * (Nheavy1_shell - Nheavy1_eff)
8518 Delta_U1 =
min(Delta_U1, 0.0);
8519 Delta_U2 = Delta_U2_shell
8520 + (Nheavy2_shell - Nheavy2_eff) * (Nheavy2_shell - Nheavy2_eff)
8522 Delta_U2 =
min(Delta_U2, 0.0);
8526 Epot0_mode1_saddle = (Nheavy1_eff - Nsymm) * (Nheavy1_eff - Nsymm) * cN_symm_sadd;
8527 Epot0_mode2_saddle = (Nheavy2_eff - Nsymm) * (Nheavy2_eff - Nsymm) * cN_symm_sadd;
8528 Epot0_symm_saddle = 0.0;
8532 Epot_mode1_saddle = Epot0_mode1_saddle + Delta_U1;
8533 Epot_mode2_saddle = Epot0_mode2_saddle + Delta_U2;
8534 Epot_symm_saddle = Epot0_symm_saddle;
8537 dUeff =
min(Epot_mode1_saddle, Epot_mode2_saddle);
8538 dUeff =
min(dUeff, Epot_symm_saddle);
8539 dUeff = dUeff - Epot_symm_saddle;
8548 epsilon0_1_saddle = Eld - Epot0_mode1_saddle;
8549 epsilon0_2_saddle = Eld - Epot0_mode2_saddle;
8552 epsilon_1_saddle = Eld - Epot_mode1_saddle;
8553 epsilon_2_saddle = Eld - Epot_mode2_saddle;
8555 epsilon_symm_saddle = Eld - Epot_symm_saddle;
8558 eexc1_saddle = epsilon_1_saddle;
8559 eexc2_saddle = epsilon_2_saddle;
8562 EEXC_MAX =
max(eexc1_saddle, eexc2_saddle);
8563 EEXC_MAX =
max(EEXC_MAX, Eld);
8566 epsilon_1_scission = Eld + E_saddle_scission - Epot_mode1_saddle;
8567 epsilon_2_scission = Eld + E_saddle_scission - Epot_mode2_saddle;
8570 epsilon_symm_scission = Eld + E_saddle_scission - Epot_symm_saddle;
8575 - Delta_U1 *
Uwash(epsilon_1_saddle /
A * Aheavy1, Ecrit, FREDSHELL, gamma_heavy1);
8577 if (E_eff1_saddle < A_levdens * hbom1 * hbom1) E_eff1_saddle = A_levdens * hbom1 * hbom1;
8579 wNasymm1_saddle = std::sqrt(
8580 0.50 * std::sqrt(1.0 / A_levdens * E_eff1_saddle)
8581 / (cN_asymm1_shell *
Uwash(epsilon_1_saddle /
A * Aheavy1, Ecrit, FREDSHELL, gamma_heavy1)
8586 - Delta_U2 *
Uwash(epsilon_2_saddle /
A * Aheavy2, Ecrit, FREDSHELL, gamma_heavy2);
8588 if (E_eff2_saddle < A_levdens * hbom2 * hbom2) E_eff2_saddle = A_levdens * hbom2 * hbom2;
8590 wNasymm2_saddle = std::sqrt(
8591 0.50 * std::sqrt(1.0 / A_levdens * E_eff2_saddle)
8592 / (cN_asymm2_shell *
Uwash(epsilon_2_saddle /
A * Aheavy2, Ecrit, FREDSHELL, gamma_heavy2)
8595 E_eff0_saddle = epsilon_symm_saddle;
8596 if (E_eff0_saddle < A_levdens * hbom3 * hbom3) E_eff0_saddle = A_levdens * hbom3 * hbom3;
8598 wNsymm_saddle = std::sqrt(0.50 * std::sqrt(1.0 / A_levdens * E_eff0_saddle) / cN_symm_sadd);
8600 if (epsilon_symm_scission > 0.0) {
8601 E_HELP =
max(E_saddle_scission, epsilon_symm_scission);
8602 wNsymm_scission = std::sqrt(0.50 * std::sqrt(1.0 / A_levdens * (E_HELP)) / cN_symm);
8605 wNsymm_scission = std::sqrt(0.50 * std::sqrt(1.0 / A_levdens * E_saddle_scission) / cN_symm);
8611 if (E_saddle_scission == 0.0) {
8612 wNasymm1_scission = wNasymm1_saddle;
8613 wNasymm2_scission = wNasymm2_saddle;
8616 if (Nheavy1_eff > 75.0) {
8617 wNasymm1_scission = std::sqrt(21.0) *
N /
A;
8618 wNasymm2_scission =
max(12.8 - 1.0 * (92.0 - Nheavy2_eff), 1.0) *
N /
A;
8621 wNasymm1_scission = wNasymm1_saddle;
8622 wNasymm2_scission = wNasymm2_saddle;
8626 wNasymm1_scission =
max(wNasymm1_scission, wNasymm1_saddle);
8627 wNasymm2_scission =
max(wNasymm2_scission, wNasymm2_saddle);
8629 wNasymm1 = wNasymm1_scission * Fwidth_asymm1;
8630 wNasymm2 = wNasymm2_scission * Fwidth_asymm2;
8631 wNsymm = wNsymm_scission * Fwidth_symm;
8634 Aheavy1_eff = Nheavy1_eff *
A /
N;
8635 Aheavy2_eff = Nheavy2_eff *
A /
N;
8637 A_levdens_heavy1 = Aheavy1_eff / xLevdens;
8638 A_levdens_heavy2 = Aheavy2_eff / xLevdens;
8639 gamma_heavy1 = A_levdens_heavy1 / (0.40 * std::pow(Aheavy1_eff, 1.3333)) * FGAMMA * FGAMMA1;
8640 gamma_heavy2 = A_levdens_heavy2 / (0.40 * std::pow(Aheavy2_eff, 1.3333)) * FGAMMA;
8642 if (epsilon_symm_saddle < A_levdens * hbom3 * hbom3)
8643 Ssymm = 2.0 * std::sqrt(A_levdens * A_levdens * hbom3 * hbom3)
8644 + (epsilon_symm_saddle - A_levdens * hbom3 * hbom3) / hbom3;
8646 Ssymm = 2.0 * std::sqrt(A_levdens * epsilon_symm_saddle);
8650 if (epsilon0_1_saddle < A_levdens * hbom1 * hbom1)
8651 Ssymm_mode1 = 2.0 * std::sqrt(A_levdens * A_levdens * hbom1 * hbom1)
8652 + (epsilon0_1_saddle - A_levdens * hbom1 * hbom1) / hbom1;
8654 Ssymm_mode1 = 2.0 * std::sqrt(A_levdens * epsilon0_1_saddle);
8656 if (epsilon0_2_saddle < A_levdens * hbom2 * hbom2)
8657 Ssymm_mode2 = 2.0 * std::sqrt(A_levdens * A_levdens * hbom2 * hbom2)
8658 + (epsilon0_2_saddle - A_levdens * hbom2 * hbom2) / hbom2;
8660 Ssymm_mode2 = 2.0 * std::sqrt(A_levdens * epsilon0_2_saddle);
8662 if (epsilon0_1_saddle
8663 - Delta_U1 *
Uwash(epsilon_1_saddle /
A * Aheavy1, Ecrit, FREDSHELL, gamma_heavy1)
8664 < A_levdens * hbom1 * hbom1)
8665 Sasymm1 = 2.0 * std::sqrt(A_levdens * A_levdens * hbom1 * hbom1)
8666 + (epsilon0_1_saddle
8667 - Delta_U1 *
Uwash(epsilon_1_saddle /
A * Aheavy1, Ecrit, FREDSHELL, gamma_heavy1)
8668 - A_levdens * hbom1 * hbom1)
8675 * (epsilon0_1_saddle
8676 - Delta_U1 *
Uwash(epsilon_1_saddle /
A * Aheavy1, Ecrit, FREDSHELL, gamma_heavy1)));
8678 if (epsilon0_2_saddle
8679 - Delta_U2 *
Uwash(epsilon_2_saddle /
A * Aheavy2, Ecrit, FREDSHELL, gamma_heavy2)
8680 < A_levdens * hbom2 * hbom2)
8681 Sasymm2 = 2.0 * std::sqrt(A_levdens * A_levdens * hbom2 * hbom2)
8682 + (epsilon0_1_saddle
8683 - Delta_U1 *
Uwash(epsilon_2_saddle /
A * Aheavy2, Ecrit, FREDSHELL, gamma_heavy2)
8684 - A_levdens * hbom2 * hbom2)
8691 * (epsilon0_2_saddle
8692 - Delta_U2 *
Uwash(epsilon_2_saddle /
A * Aheavy2, Ecrit, FREDSHELL, gamma_heavy2)));
8694 Yasymm1 = (std::exp(Sasymm1 - Ssymm) - std::exp(Ssymm_mode1 - Ssymm)) * wNasymm1_saddle
8695 / wNsymm_saddle * 2.0;
8697 Yasymm2 = (std::exp(Sasymm2 - Ssymm) - std::exp(Ssymm_mode2 - Ssymm)) * wNasymm2_saddle
8698 / wNsymm_saddle * 2.0;
8700 Ysum = Ysymm + Yasymm1 + Yasymm2;
8703 Ysymm = Ysymm / Ysum;
8704 Yasymm1 = Yasymm1 / Ysum;
8705 Yasymm2 = Yasymm2 / Ysum;
8706 Yasymm = Yasymm1 + Yasymm2;
8713 if ((epsilon_symm_saddle < epsilon_1_saddle) && (epsilon_symm_saddle < epsilon_2_saddle))
8715 else if (epsilon_1_saddle < epsilon_2_saddle)
8722 if (
mod(Z, 2.0) == 0)
8723 r_e_o = std::pow(10.0, -0.0170 * (E_saddle_scission + Eld) * (E_saddle_scission + Eld));
8741 if (rmode < Yasymm1)
8743 else if ((rmode > Yasymm1) && (rmode < Yasymm))
8752 N1mean = Nheavy1_eff;
8757 N1mean = Nheavy2_eff;
8774 while (N1r < 5.0 || N2r < 5.0) {
8777 N1r =
gausshaz(0, N1mean, N1width);
8784 Z1UCD = Z /
N * N1r;
8785 Z2UCD = Z /
N * N2r;
8793 E_scission_pre =
max(epsilon_1_scission, 1.0);
8796 if (N1mean >
N * 0.50) {
8807 E_scission_pre =
max(epsilon_2_scission, 1.0);
8808 if (N1mean >
N * 0.50) {
8809 beta1 = (N1r - 92.0) * 0.030 + 0.60;
8811 beta1gs = ecld->beta2[
idint(N1r)][
idint(Z1UCD)];
8812 beta2gs = ecld->beta2[
idint(N2r)][
idint(Z2UCD)];
8814 beta1 =
max(beta1, beta1gs);
8815 beta2 = 1.0 - beta1;
8816 beta2 =
max(beta2, beta2gs);
8819 beta1gs = ecld->beta2[
idint(N1r)][
idint(Z1UCD)];
8820 beta2gs = ecld->beta2[
idint(N2r)][
idint(Z2UCD)];
8822 beta2 = (N2r - 92.0) * 0.030 + 0.60;
8823 beta2 =
max(beta2, beta2gs);
8824 beta1 = 1.0 - beta2;
8825 beta1 =
max(beta1, beta1gs);
8837 betags = ecld->beta2[
idint(Nsymm)][
idint(Zsymm)];
8838 beta1gs = ecld->beta2[
idint(N1r)][
idint(Z1UCD)];
8839 beta2gs = ecld->beta2[
idint(N2r)][
idint(Z2UCD)];
8840 beta =
max(0.177963 + 0.0153241 * Zsymm - 1.62037e-4 * Zsymm * Zsymm, betags);
8841 beta1 =
max(0.177963 + 0.0153241 * Z1UCD - 1.62037e-4 * Z1UCD * Z1UCD, beta1gs);
8842 beta2 =
max(0.177963 + 0.0153241 * Z2UCD - 1.62037e-4 * Z2UCD * Z2UCD, beta2gs);
8844 E_asym =
frldm(Z1UCD, N1r, beta1) +
frldm(Z2UCD, N2r, beta2)
8845 +
ecoul(Z1UCD, N1r, beta1, Z2UCD, N2r, beta2, 2.0) - 2.0 *
frldm(Zsymm, Nsymm, beta)
8846 -
ecoul(Zsymm, Nsymm, beta, Zsymm, Nsymm, beta, 2.0);
8847 E_scission_pre =
max(epsilon_symm_scission - E_asym, 1.);
8856 if (E_scission_pre > 5. && NbLam0 < 1) {
8857 evap_postsaddle(
A, Z, E_scission_pre, &E_scission_post, &A_scission, &Z_scission, vx_eva_sc,
8858 vy_eva_sc, vz_eva_sc, &NbLam0);
8859 N_scission = A_scission - Z_scission;
8864 E_scission_post = E_scission_pre;
8865 N_scission = A_scission - Z_scission;
8871 N1r = N1r * N_scission /
N;
8872 N2r = N2r * N_scission /
N;
8873 Z1UCD = Z1UCD * Z_scission / Z;
8874 Z2UCD = Z2UCD * Z_scission / Z;
8911 CZ = (
frldm(Z1UCD - 1.0, N1r + 1.0, beta1) +
frldm(Z2UCD + 1.0, N2r - 1.0, beta2)
8912 +
frldm(Z1UCD + 1.0, N1r - 1.0, beta1) +
frldm(Z2UCD - 1.0, N2r + 1.0, beta2)
8913 +
ecoul(Z1UCD - 1.0, N1r + 1.0, beta1, Z2UCD + 1.0, N2r - 1.0, beta2, 2.0)
8914 +
ecoul(Z1UCD + 1.0, N1r - 1.0, beta1, Z2UCD - 1.0, N2r + 1.0, beta2, 2.0)
8915 - 2.0 *
ecoul(Z1UCD, N1r, beta1, Z2UCD, N2r, beta2, 2.0) - 2.0 *
frldm(Z1UCD, N1r, beta1)
8916 - 2.0 *
frldm(Z2UCD, N2r, beta2))
8919 if (1.0 / A_levdens * E_scission_post < 0.0)
8920 std::cout <<
"DSQRT 1 < 0" << A_levdens <<
" " << E_scission_post << std::endl;
8922 if (0.50 * std::sqrt(1.0 / A_levdens * E_scission_post) / CZ < 0.0) {
8923 std::cout <<
"DSQRT 2 < 0 " << CZ << std::endl;
8924 std::cout <<
"This event was not considered" << std::endl;
8928 ZA1width = std::sqrt(0.5 * std::sqrt(1.0 / A_levdens * E_scission_post) / CZ);
8939 ZA1width =
max(ZA1width, sigZmin);
8941 if (imode == 1 && cpol1 != 0.0) {
8945 Z1rr = Z1UCD - cpol1 * A_scission / N_scission;
8951 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN "
8952 "CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED"
8956 if ((
utilabs(Z1rr - Z1r) > 3.0 * ZA1width) || Z1r < 1.0)
goto fiss2801;
8960 if (imode == 2 && cpol2 != 0.0) {
8964 Z1rr = Z1UCD - cpol2 * A_scission / N_scission;
8969 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN "
8970 "CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED"
8974 if ((
utilabs(Z1rr - Z1r) > 3.0 * ZA1width) || Z1r < 1.0)
goto fiss2802;
8984 re1 =
frldm(Z1UCD - 1.0, N1r + 1.0, beta1) +
frldm(Z2UCD + 1.0, N2r - 1.0, beta2)
8985 +
ecoul(Z1UCD - 1.0, N1r + 1.0, beta1, Z2UCD + 1.0, N2r - 1.0, beta2, d);
8986 re2 =
frldm(Z1UCD, N1r, beta1) +
frldm(Z2UCD, N2r, beta2)
8987 +
ecoul(Z1UCD, N1r, beta1, Z2UCD, N2r, beta2, d);
8988 re3 =
frldm(Z1UCD + 1.0, N1r - 1.0, beta1) +
frldm(Z2UCD - 1.0, N2r + 1.0, beta2)
8989 +
ecoul(Z1UCD + 1.0, N1r - 1.0, beta1, Z2UCD - 1.0, N2r + 1.0, beta2, d);
8990 eps2 = (re1 - 2.0 * re2 + re3) / 2.0;
8991 eps1 = (re3 - re1) / 2.0;
8992 DN1_POL = -eps1 / (2.0 * eps2);
8994 Z1rr = Z1UCD + DN1_POL;
8999 DN1_POL = DN1_POL - 0.6 *
Uwash(E_scission_post, Ecrit, FREDSHELL, gamma);
9000 Z1rr = Z1UCD + DN1_POL;
9001 if (Z1rr < 50.) Z1rr = 50.0;
9004 DN1_POL = DN1_POL + 0.60 *
Uwash(E_scission_post, Ecrit, FREDSHELL, gamma);
9005 Z1rr = Z1UCD + DN1_POL;
9006 if (Z1rr > 50.0) Z1rr = 50.0;
9016 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN "
9017 "CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED"
9022 if ((
utilabs(Z1rr - Z1r) > 3.0 * ZA1width) || (Z1r < 1.0))
goto fiss2803;
9033 z2 =
dint(Z_scission) - z1;
9035 N2 =
dint(N_scission) - N1;
9039 if ((z1 < 0) || (z2 < 0) || (a1 < 0) || (a2 < 0)) {
9040 std::cout <<
" -------------------------------" << std::endl;
9041 std::cout <<
" Z, A, N : " << Z <<
" " <<
A <<
" " <<
N << std::endl;
9042 std::cout << z1 <<
" " << z2 <<
" " << a1 <<
" " << a2 << std::endl;
9043 std::cout << E_scission_post <<
" " << A_levdens <<
" " << CZ << std::endl;
9045 std::cout <<
" -------------------------------" << std::endl;
9054 if (N1mean >
N * 0.50) {
9058 if (beta2 < beta2gs) beta2 = beta2gs;
9059 E1exc = E_scission_pre * a1 /
A + E_defo;
9060 E_defo =
frldm(z2, N2, beta2) -
frldm(z2, N2, beta2gs);
9061 E2exc = E_scission_pre * a2 /
A + E_defo;
9066 if (beta1 < beta1gs) beta1 = beta1gs;
9067 E_defo =
frldm(z1, N1, beta1) -
frldm(z1, N1, beta1gs);
9068 E1exc = E_scission_pre * a1 /
A + E_defo;
9070 E2exc = E_scission_pre * a2 /
A + E_defo;
9076 if (N1mean >
N * 0.5) {
9079 if (beta1 < beta1gs) beta1 = beta1gs;
9080 E_defo =
frldm(z1, N1, beta1) -
frldm(z1, N1, beta1gs);
9081 E1exc = E_scission_pre * a1 /
A + E_defo;
9083 if (beta2 < beta2gs) beta2 = beta2gs;
9084 E_defo =
frldm(z2, N2, beta2) -
frldm(z2, N2, beta2gs);
9085 E2exc = E_scission_pre * a2 /
A + E_defo;
9090 if (beta2 < beta2gs) beta2 = beta2gs;
9091 E_defo =
frldm(z2, N2, beta2) -
frldm(z2, N2, beta2gs);
9092 E2exc = E_scission_pre * a2 /
A + E_defo;
9094 if (beta1 < beta1gs) beta1 = beta1gs;
9095 E_defo =
frldm(z1, N1, beta1) -
frldm(z1, N1, beta1gs);
9096 E1exc = E_scission_pre * a1 /
A + E_defo;
9103 if (beta1 < beta1gs) beta1 = beta1gs;
9105 if (beta2 < beta2gs) beta2 = beta2gs;
9106 E_defo1 =
frldm(z1, N1, beta1) -
frldm(z1, N1, beta1gs);
9107 E_defo2 =
frldm(z2, N2, beta2) -
frldm(z2, N2, beta2gs);
9108 E1exc = E_scission_pre * a1 /
A + E_defo1;
9109 E2exc = E_scission_pre * a2 /
A + E_defo2;
9113 TKER = (z1 * z2 * 1.440)
9114 / (R0 * std::pow(a1, 0.333330) * (1.0 + 2.0 / 3.0 * beta1)
9115 + R0 * std::pow(a2, 0.333330) * (1.0 + 2.0 / 3.0 * beta2) + 2.0);
9117 EkinR1 = TKER * a2 /
A;
9118 EkinR2 = TKER * a1 /
A;
9119 v1 = std::sqrt(EkinR1 / a1) * 1.3887;
9120 v2 = std::sqrt(EkinR2 / a2) * 1.3887;
9128 e1 =
gausshaz(0, E1exc, E1exc_sigma);
9129 if (e1 < 0.)
goto fis987;
9132 e2 =
gausshaz(0, E2exc, E2exc_sigma);
9133 if (e2 < 0.)
goto fis988;
9135 (*NbLam0_par) = NbLam0;
9420 G4int INMIN, INMAX, NDIF = 0, IMEM;
9421 G4int NEVA = 0, PEVA = 0;
9430 BU_TAB_TEMP[i][0] = 0.0;
9431 BU_TAB_TEMP[i][1] = 0.0;
9432 BU_TAB_TEMP[i][2] = 0.0;
9433 BU_TAB_TEMP[i][3] = 0.0;
9434 BU_TAB_TEMP[i][4] = 0.0;
9441 if (AFP == 0 && ZFP == 0) {
9445 if ((AFP == 1 && ZFP == 0) || (AFP == 1 && ZFP == 1) || (AFP == 2 && ZFP == 1)
9446 || (AFP == 3 && ZFP == 1) || (AFP == 3 && ZFP == 2) || (AFP == 4 && ZFP == 2)
9447 || (AFP == 6 && ZFP == 2) || (AFP == 8 && ZFP == 2))
9455 if ((AFP - ZFP) == 0 && ZFP > 1) {
9456 for (
G4int I = 0; I <= AFP - 2; I++) {
9458 G4double(AFP - I - 1), VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y,
9460 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9461 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9462 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9463 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9464 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9465 *ILOOP = *ILOOP + 1;
9482 for (
G4int I = 1; I <= 10; I++) {
9484 if (INMIN <= NDIF) {
9494 for (
G4int I = 0; I < IMEM; I++) {
9496 G4double(NDIF + ZFP + IMEM - I - 1),
G4double(ZFP + IMEM - I - 1), VX, VY, VZ,
9497 &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9498 BU_TAB_TEMP[I + 1 + *ILOOP][0] = 1.0;
9499 BU_TAB_TEMP[I + 1 + *ILOOP][1] = 1.0;
9500 BU_TAB_TEMP[I + 1 + *ILOOP][2] = VP2X;
9501 BU_TAB_TEMP[I + 1 + *ILOOP][3] = VP2Y;
9502 BU_TAB_TEMP[I + 1 + *ILOOP][4] = VP2Z;
9507 *ILOOP = *ILOOP + IMEM;
9511 NEVA = NDIF - INMAX;
9514 for (
G4int I = 0; I < NEVA; I++) {
9516 G4double(ZFP), VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9518 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9519 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9520 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9521 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9522 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9523 *ILOOP = *ILOOP + 1;
9530 if ((AFP >= 2) && (ZFP == 0)) {
9531 for (
G4int I = 0; I <= AFP - 2; I++) {
9533 VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9535 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9536 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9537 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9538 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9539 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9540 *ILOOP = *ILOOP + 1;
9552 std::cout <<
"WARNING - BU UNSTABLE: AF < ZF" << std::endl;
9557 if ((AFP >= 4) && (ZFP == 1)) {
9560 for (
G4int I = 0; I < AFP - 3; I++) {
9562 VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9564 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9565 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9566 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9567 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9568 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9569 *ILOOP = *ILOOP + 1;
9581 if ((AFP == 4) && (ZFP == 3)) {
9587 unstable_tke(4.0, 3.0, 3.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9589 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9590 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9591 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9592 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9593 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9594 *ILOOP = *ILOOP + 1;
9596 if ((AFP == 5) && (ZFP == 2)) {
9602 unstable_tke(5.0, 2.0, 4.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9603 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9604 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9605 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9606 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9607 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9608 *ILOOP = *ILOOP + 1;
9611 if ((AFP == 5) && (ZFP == 3)) {
9617 unstable_tke(5.0, 3.0, 4.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9618 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9619 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9620 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9621 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9622 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9623 *ILOOP = *ILOOP + 1;
9626 if ((AFP == 6) && (ZFP == 4)) {
9633 unstable_tke(6.0, 4.0, 5.0, 3.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9634 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9635 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9636 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9637 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9638 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9639 *ILOOP = *ILOOP + 1;
9645 unstable_tke(5.0, 3.0, 4.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9646 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9647 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9648 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9649 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9650 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9651 *ILOOP = *ILOOP + 1;
9653 if ((AFP == 7) && (ZFP == 2)) {
9659 unstable_tke(7.0, 2.0, 6.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9660 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9661 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9662 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9663 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9664 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9665 *ILOOP = *ILOOP + 1;
9668 if ((AFP == 7) && (ZFP == 5)) {
9670 for (
int I = 0; I <= AFP - 5; I++) {
9671 unstable_tke(
double(AFP - I),
double(ZFP - I),
double(AFP - I - 1),
double(ZFP - I - 1), VX,
9672 VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9673 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9674 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9675 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9676 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9677 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9678 *ILOOP = *ILOOP + 1;
9689 if ((AFP == 8) && (ZFP == 4)) {
9694 unstable_tke(8.0, 4.0, 4.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9695 BU_TAB_TEMP[*ILOOP][0] = 2.0;
9696 BU_TAB_TEMP[*ILOOP][1] = 4.0;
9697 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9698 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9699 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9700 *ILOOP = *ILOOP + 1;
9702 if ((AFP == 8) && (ZFP == 6)) {
9709 unstable_tke(8.0, 6.0, 7.0, 5.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9710 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9711 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9712 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9713 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9714 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9715 *ILOOP = *ILOOP + 1;
9720 unstable_tke(7.0, 5.0, 6.0, 4.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9721 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9722 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9723 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9724 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9725 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9726 *ILOOP = *ILOOP + 1;
9732 if ((AFP == 9) && (ZFP == 2)) {
9739 unstable_tke(9.0, 2.0, 8.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9740 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9741 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9742 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9743 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9744 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9745 *ILOOP = *ILOOP + 1;
9751 if ((AFP == 9) && (ZFP == 5)) {
9757 unstable_tke(9.0, 5.0, 8.0, 4.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9758 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9759 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9760 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9761 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9762 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9763 *ILOOP = *ILOOP + 1;
9768 unstable_tke(8.0, 4.0, 4.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9769 BU_TAB_TEMP[*ILOOP][0] = 2.0;
9770 BU_TAB_TEMP[*ILOOP][1] = 4.0;
9771 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9772 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9773 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9774 *ILOOP = *ILOOP + 1;
9780 if ((AFP == 10) && (ZFP == 2)) {
9787 unstable_tke(10.0, 2.0, 9.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y,
9789 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9790 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9791 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9792 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9793 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9794 *ILOOP = *ILOOP + 1;
9800 unstable_tke(9.0, 2.0, 8.0, 2.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
9801 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9802 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9803 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9804 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9805 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9806 *ILOOP = *ILOOP + 1;
9811 if ((AFP == 10) && (ZFP == 3)) {
9817 unstable_tke(10.0, 3.0, 9.0, 3.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y,
9819 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9820 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9821 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9822 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9823 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9824 *ILOOP = *ILOOP + 1;
9829 if ((AFP == 10) && (ZFP == 7)) {
9835 unstable_tke(10.0, 7.0, 9.0, 6.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y,
9837 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9838 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9839 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9840 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9841 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9842 *ILOOP = *ILOOP + 1;
9848 if ((AFP == 11) && (ZFP == 7)) {
9854 unstable_tke(11.0, 7.0, 10.0, 6.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y,
9856 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9857 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9858 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9859 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9860 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9861 *ILOOP = *ILOOP + 1;
9866 if ((AFP == 12) && (ZFP == 8)) {
9873 unstable_tke(12.0, 8.0, 11.0, 7.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y,
9875 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9876 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9877 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9878 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9879 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9880 *ILOOP = *ILOOP + 1;
9885 unstable_tke(11.0, 7.0, 10.0, 6.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y,
9887 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9888 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9889 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9890 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9891 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9892 *ILOOP = *ILOOP + 1;
9897 if ((AFP == 15) && (ZFP == 9)) {
9903 unstable_tke(15.0, 9.0, 14.0, 8.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y,
9905 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9906 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9907 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9908 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9909 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9910 *ILOOP = *ILOOP + 1;
9916 if ((AFP == 16) && (ZFP == 9)) {
9922 unstable_tke(16.0, 9.0, 15.0, 8.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y,
9924 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9925 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9926 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9927 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9928 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9929 *ILOOP = *ILOOP + 1;
9935 if ((AFP == 16) && (ZFP == 10)) {
9941 unstable_tke(16.0, 10.0, 15.0, 9.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y,
9943 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9944 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9945 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9946 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9947 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9948 *ILOOP = *ILOOP + 1;
9953 unstable_tke(15.0, 9.0, 14.0, 8.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y,
9955 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9956 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9957 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9958 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9959 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9960 *ILOOP = *ILOOP + 1;
9965 if ((AFP == 18) && (ZFP == 11)) {
9971 unstable_tke(18.0, 11.0, 17.0, 10.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y,
9973 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9974 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9975 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9976 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9977 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9978 *ILOOP = *ILOOP + 1;
9983 if ((AFP == 19) && (ZFP == 11)) {
9989 unstable_tke(19.0, 11.0, 18.0, 10.0, VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y,
9991 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9992 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9993 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9994 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9995 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9996 *ILOOP = *ILOOP + 1;
10001 if (ZFP >= 4 && (AFP - ZFP) == 1) {
10006 for (
G4int I = 0; I < NEVA; I++) {
10008 VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y, &VP2Z);
10009 BU_TAB_TEMP[*ILOOP][0] = 0.0;
10010 BU_TAB_TEMP[*ILOOP][1] = 1.0;
10011 BU_TAB_TEMP[*ILOOP][2] = VP2X;
10012 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
10013 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
10014 *ILOOP = *ILOOP + 1;
10019 for (
G4int I = 0; I < PEVA; I++) {
10021 G4double(ZFP - I - 1), VX, VY, VZ, &(*VP1X), &(*VP1Y), &(*VP1Z), &VP2X, &VP2Y,
10023 BU_TAB_TEMP[*ILOOP][0] = 1.0;
10024 BU_TAB_TEMP[*ILOOP][1] = 1.0;
10025 BU_TAB_TEMP[*ILOOP][2] = VP2X;
10026 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
10027 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
10028 *ILOOP = *ILOOP + 1;