32 #define ABLAXX_IN_GEANT4_MODE 1
52 fissionDistri(A,Z,E,A1,Z1,E1,K1,A2,Z2,E2,K2);
88 G4double r_in = 0.0, r_rest = 0.0, r_help = 0.0;
94 r_in = r_origin + 0.5;
95 r_floor = (float)((
int)(r_in));
96 if (r_even_odd < 0.001) {
97 i_out = (int)(r_floor);
100 r_rest = r_in - r_floor;
101 r_middle = r_floor + 0.5;
102 n_floor = (int)(r_floor);
103 if (n_floor%2 == 0) {
105 r_help = r_middle + (r_rest - 0.5) * (1.0 - r_even_odd);
109 r_help = r_middle + (r_rest - 0.5) * (1.0 + r_even_odd);
111 i_out = (int)(r_help);
127 G4double xcom = 0.0, xvs = 0.0, xe = 0.0;
131 alpha = ( std::sqrt(5.0/(4.0*pi)) ) * beta;
133 xcom = 1.0 - 1.7826 * ((a - 2.0*
z)/a)*((a - 2.0*
z)/a);
135 xvs = - xcom * ( 15.4941 * a -
136 17.9439 * std::pow(a,2.0/3.0) * (1.0+0.4*alpha*
alpha) );
138 xe = z*z * (0.7053/(std::pow(a,1.0/3.0)) * (1.0-0.2*alpha*alpha) - 1.1529/
a);
164 dtot = r0 * ( std::pow((z1+n1),1.0/3.0) * (1.0+0.6666*beta1)
165 + std::pow((z2+n2),1.0/3.0) * (1.0+0.6666*beta2) ) + d;
166 fecoul = z1 * z2 * 1.44 / dtot;
249 G4double aheavy1 = 0.0, aheavy2 = 0.0;
251 G4double zheavy1_shell = 0.0, zheavy2_shell = 0.0;
253 G4double sasymm1 = 0.0, sasymm2 = 0.0, ssymm = 0.0, ysum = 0.0;
254 G4double ssymm_mode1 = 0.0, ssymm_mode2 = 0.0;
256 G4double wzasymm1_saddle, wzasymm2_saddle, wzsymm_saddle = 0.0;
257 G4double wzasymm1_scission = 0.0, wzasymm2_scission = 0.0, wzsymm_scission = 0.0;
258 G4double wzasymm1 = 0.0, wzasymm2 = 0.0, wzsymm = 0.0;
261 G4double z1mean = 0.0, z1width = 0.0;
265 G4double zsymm = 0.0, nsymm = 0.0, asymm = 0.0;
266 G4double n1mean = 0.0, n1width = 0.0;
269 G4double re1 = 0.0, re2 = 0.0, re3 = 0.0;
279 G4double a_levdens_heavy1 = 0.0, a_levdens_heavy2 = 0.0;
282 G4double epsilon_1_saddle = 0.0, epsilon0_1_saddle = 0.0;
283 G4double epsilon_2_saddle = 0.0, epsilon0_2_saddle = 0.0, epsilon_symm_saddle = 0.0;
284 G4double epsilon_1_scission = 0.0, epsilon0_1_scission = 0.0;
285 G4double epsilon_2_scission = 0.0, epsilon0_2_scission = 0.0;
286 G4double epsilon_symm_scission = 0.0;
288 G4double e_eff1_saddle = 0.0, e_eff2_saddle = 0.0;
289 G4int icz = 0, k = 0;
293 G4double ed1_low = 0.0, ed2_low = 0.0, ed1_high = 0.0, ed2_high = 0.0, ed1 = 0.0, ed2 = 0.0;
316 const G4double delta_u1_shell = -2.5;
319 const G4double delta_u2_shell = -5.5;
322 const G4double cz_asymm1_shell = 0.7;
324 const G4double cz_asymm2_shell = 0.08;
344 const G4double friction_factor = 1.0;
376 const G4int itest = 0;
378 G4double reps1 = 0.0, reps2 = 0.0, rn1_pol = 0.0;
387 const G4int maxCount = 500;
394 if ( (std::pow(z,2)/a < 25.0) || (n < nheavy2) || (e > 500.0) ) {
406 zheavy1_shell = ((nheavy1/
n) * z) - 0.8;
408 zheavy2_shell = ((nheavy2/
n) * z) - 0.8;
414 e_saddle_scission = (3.535 * std::pow(z,2)/a - 121.1) * friction_factor;
420 if (e_saddle_scission < 0.0) {
421 e_saddle_scission = 0.0;
435 if ( (std::pow(z,2))/a < 33.9186) {
436 masscurv = std::pow( 10.0,(-1.093364 + 0.082933 * (std::pow(z,2)/a)
437 - 0.0002602 * (std::pow(z,4)/std::pow(a,2))) );
439 masscurv = std::pow( 10.0,(3.053536 - 0.056477 * (std::pow(z,2)/a)
440 + 0.0002454 * (std::pow(z,4)/std::pow(a,2))) );
443 cz_symm = (8.0/std::pow(z,2)) * masscurv;
459 asymm = nsymm + zsymm;
461 zheavy1 = (cz_symm*zsymm + cz_asymm1_shell*zheavy1_shell)/(cz_symm + cz_asymm1_shell);
462 zheavy2 = (cz_symm*zsymm + cz_asymm2_shell*zheavy2_shell)/(cz_symm + cz_asymm2_shell);
465 nheavy1_eff = (zheavy1 + 0.8)*(n/
z);
466 nheavy2_eff = (zheavy2 + 0.8)*(n/
z);
467 aheavy1 = nheavy1_eff + zheavy1;
468 aheavy2 = nheavy2_eff + zheavy2;
474 a_levdens_heavy1 = aheavy1 / 8.0;
475 a_levdens_heavy2 = aheavy2 / 8.0;
476 gamma = a_levdens / (0.4 * (std::pow(a,1.3333)) );
477 gamma_heavy1 = ( a_levdens_heavy1 / (0.4 * (std::pow(aheavy1,1.3333)) ) ) * fgamma1;
478 gamma_heavy2 = a_levdens_heavy2 / (0.4 * (std::pow(aheavy2,1.3333)) );
482 cn =
umass(zsymm,(nsymm+1.),0.0) +
umass(zsymm,(nsymm-1.),0.0)
483 + 1.44 * (std::pow(zsymm,2))/
484 ( (std::pow(r_null,2)) *
485 ( std::pow((asymm+1.0),1.0/3.0) + std::pow((asymm-1.0),1.0/3.0) ) *
486 ( std::pow((asymm+1.0),1.0/3.0) + std::pow((asymm-1.0),1.0/3.0) ) )
487 - 2.0 *
umass(zsymm,nsymm,0.0)
488 - 1.44 * (std::pow(zsymm,2))/
489 ( ( 2.0 * r_null * (std::pow(asymm,1.0/3.0)) ) *
490 ( 2.0 * r_null * (std::pow(asymm,1.0/3.0)) ) );
493 delta_u1 = delta_u1_shell + (std::pow((zheavy1_shell-zheavy1),2))*cz_asymm1_shell;
495 delta_u2 = delta_u2_shell + (std::pow((zheavy2_shell-zheavy2),2))*cz_asymm2_shell;
498 Basym_1 = Bsym + std::pow((zheavy1-zsymm), 2) * cz_symm + delta_u1;
499 Basym_2 = Bsym + std::pow((zheavy2-zsymm), 2) * cz_symm + delta_u2;
500 if(Bsym < Basym_1 && Bsym < Basym_2) {
503 epsilon0_1_saddle = (e + e_zero_point - std::pow((zheavy1 - zsymm), 2) * cz_symm);
504 epsilon0_2_saddle = (e + e_zero_point - std::pow((zheavy2 - zsymm), 2) * cz_symm);
506 epsilon_1_saddle = epsilon0_1_saddle - delta_u1;
507 epsilon_2_saddle = epsilon0_2_saddle - delta_u2;
509 epsilon_symm_saddle = e + e_zero_point;
510 eexc1_saddle = epsilon_1_saddle;
511 eexc2_saddle = epsilon_2_saddle;
515 epsilon0_1_scission = (e + e_saddle_scission - std::pow((zheavy1 - zsymm), 2) * cz_symm) * aheavy1/a;
516 epsilon_1_scission = epsilon0_1_scission - delta_u1*(aheavy1/
a);
518 epsilon0_2_scission = (e + e_saddle_scission - std::pow((zheavy2 - zsymm), 2) * cz_symm) * aheavy2/a;
519 epsilon_2_scission = epsilon0_2_scission - delta_u2*(aheavy2/
a);
521 epsilon_symm_scission = e + e_saddle_scission;
522 }
else if (Basym_1 < Bsym && Basym_1 < Basym_2) {
525 epsilon_symm_saddle = (e + e_zero_point + delta_u1 + std::pow((zheavy1-zsymm), 2) * cz_symm);
526 epsilon0_2_saddle = (epsilon_symm_saddle - std::pow((zheavy2-zsymm), 2) * cz_symm);
527 epsilon_2_saddle = epsilon0_2_saddle - delta_u2;
528 epsilon0_1_saddle = e + e_zero_point + delta_u1;
529 epsilon_1_saddle = e + e_zero_point;
530 eexc1_saddle = epsilon_1_saddle;
531 eexc2_saddle = epsilon_2_saddle;
535 epsilon_symm_scission = (e + e_saddle_scission + std::pow((zheavy1-zsymm), 2) * cz_symm + delta_u1);
536 epsilon0_2_scission = (epsilon_symm_scission - std::pow((zheavy2-zsymm), 2) * cz_symm) * aheavy2/a;
537 epsilon_2_scission = epsilon0_2_scission - delta_u2*aheavy2/
a;
538 epsilon0_1_scission = (e + e_saddle_scission + delta_u1) * aheavy1/a;
539 epsilon_1_scission = (e + e_saddle_scission) * aheavy1/a;
540 }
else if (Basym_2 < Bsym && Basym_2 < Basym_1) {
543 epsilon_symm_saddle = (e + e_zero_point + delta_u2 + std::pow((zheavy2-zsymm), 2) * cz_symm);
544 epsilon0_1_saddle = (epsilon_symm_saddle - std::pow((zheavy1-zsymm), 2) * cz_symm);
545 epsilon_1_saddle = epsilon0_1_saddle - delta_u1;
546 epsilon0_2_saddle = e + e_zero_point + delta_u2;
547 epsilon_2_saddle = e + e_zero_point;
548 eexc1_saddle = epsilon_1_saddle;
549 eexc2_saddle = epsilon_2_saddle;
553 epsilon_symm_scission = (e + e_saddle_scission + std::pow((zheavy2-zsymm), 2) * cz_symm + delta_u2);
554 epsilon0_1_scission = (epsilon_symm_scission - std::pow((zheavy1-zsymm), 2) * cz_symm) * aheavy1/a;
555 epsilon_1_scission = epsilon0_1_scission - delta_u1*aheavy1/
a;
556 epsilon0_2_scission = (e + e_saddle_scission + delta_u2) * aheavy2/a;
557 epsilon_2_scission = (e + e_saddle_scission) * aheavy2/a;
562 if(epsilon_1_saddle < 0.0) epsilon_1_saddle = 0.0;
563 if(epsilon_2_saddle < 0.0) epsilon_2_saddle = 0.0;
564 if(epsilon0_1_saddle < 0.0) epsilon0_1_saddle = 0.0;
565 if(epsilon0_2_saddle < 0.0) epsilon0_2_saddle = 0.0;
566 if(epsilon_symm_saddle < 0.0) epsilon_symm_saddle = 0.0;
568 if(epsilon_1_scission < 0.0) epsilon_1_scission = 0.0;
569 if(epsilon_2_scission < 0.0) epsilon_2_scission = 0.0;
570 if(epsilon0_1_scission < 0.0) epsilon0_1_scission = 0.0;
571 if(epsilon0_2_scission < 0.0) epsilon0_2_scission = 0.0;
572 if(epsilon_symm_scission < 0.0) epsilon_symm_scission = 0.0;
578 e_eff1_saddle = epsilon0_1_saddle - delta_u1 * (std::exp((-epsilon_1_saddle*gamma)));
580 if (e_eff1_saddle > 0.0) {
581 wzasymm1_saddle = std::sqrt( (0.5) *
582 (std::sqrt(1.0/a_levdens*e_eff1_saddle)) /
583 (cz_asymm1_shell * std::exp((-epsilon_1_saddle*gamma)) + cz_symm) );
585 wzasymm1_saddle = 1.0;
588 e_eff2_saddle = epsilon0_2_saddle - delta_u2 * std::exp((-epsilon_2_saddle*gamma));
589 if (e_eff2_saddle > 0.0) {
590 wzasymm2_saddle = std::sqrt( (0.5 *
591 (std::sqrt(1.0/a_levdens*e_eff2_saddle)) /
592 (cz_asymm2_shell * std::exp((-epsilon_2_saddle*gamma)) + cz_symm) ) );
594 wzasymm2_saddle = 1.0;
597 if(e - e_zero_point > 0.0) {
598 wzsymm_saddle = std::sqrt( (0.5 *
599 (std::sqrt(1.0/a_levdens*(e+epsilon_symm_saddle))) / cz_symm ) );
613 wzsymm_scission = wzsymm_saddle;
615 if (e_saddle_scission == 0.0) {
616 wzasymm1_scission = wzasymm1_saddle;
617 wzasymm2_scission = wzasymm2_saddle;
619 if (nheavy1_eff > 75.0) {
620 wzasymm1_scission = (std::sqrt(21.0)) * z/a;
621 double RR = (70.0-28.0)/3.0*(z*z/a-35.0)+28.0;
623 wzasymm2_scission = std::sqrt(RR)*(z/
a);
625 wzasymm2_scission = 0.0;
628 wzasymm1_scission = wzasymm1_saddle;
629 wzasymm2_scission = wzasymm2_saddle;
633 wzasymm1_scission =
max(wzasymm1_scission,wzasymm1_saddle);
634 wzasymm2_scission =
max(wzasymm2_scission,wzasymm2_saddle);
636 wzasymm1 = wzasymm1_scission * fwidth_asymm1;
637 wzasymm2 = wzasymm2_scission * fwidth_asymm2;
638 wzsymm = wzsymm_scission;
663 n1width = std::sqrt(0.5 * std::sqrt(1.0/a_levdens*(e + e_saddle_scission)) / cn + 1.35);
664 if ( (epsilon0_1_saddle - delta_u1*std::exp((-epsilon_1_saddle*gamma_heavy1))) < 0.0) {
667 sasymm1 = 2.0 * std::sqrt( a_levdens * (epsilon0_1_saddle -
668 delta_u1*(std::exp((-epsilon_1_saddle*gamma_heavy1))) ) );
671 if ( (epsilon0_2_saddle - delta_u2*std::exp((-epsilon_2_saddle*gamma_heavy2))) < 0.0) {
674 sasymm2 = 2.0 * std::sqrt( a_levdens * (epsilon0_2_saddle -
675 delta_u2*(std::exp((-epsilon_2_saddle*gamma_heavy2))) ) );
678 if (epsilon_symm_saddle > 0.0) {
679 ssymm = 2.0 * std::sqrt( a_levdens*(epsilon_symm_saddle) );
686 if (epsilon0_1_saddle < 0.0) {
687 yasymm1 = std::exp((sasymm1-ssymm)) * wzasymm1_saddle / wzsymm_saddle * 2.0;
690 ssymm_mode1 = 2.0 * std::sqrt( a_levdens*(epsilon0_1_saddle) );
691 yasymm1 = ( std::exp((sasymm1-ssymm)) - std::exp((ssymm_mode1 - ssymm)) )
692 * wzasymm1_saddle / wzsymm_saddle * 2.0;
695 if (epsilon0_2_saddle < 0.0) {
696 yasymm2 = std::exp((sasymm2-ssymm)) * wzasymm2_saddle / wzsymm_saddle * 2.0;
699 ssymm_mode2 = 2.0 * std::sqrt( a_levdens*(epsilon0_2_saddle) );
700 yasymm2 = ( std::exp((sasymm2-ssymm)) - std::exp((ssymm_mode2 - ssymm)) )
701 * wzasymm2_saddle / wzsymm_saddle * 2.0;
707 if ( (sasymm1 > -10.0) && (sasymm2 > -10.0) ) {
709 yasymm1 = std::exp(sasymm1) * wzasymm1_saddle * 2.0;
710 yasymm2 = std::exp(sasymm2) * wzasymm2_saddle * 2.0;
715 ysum = ysymm + yasymm1 + yasymm2;
717 ysymm = ysymm / ysum;
718 yasymm1 = yasymm1 / ysum;
719 yasymm2 = yasymm2 / ysum;
725 if ( (epsilon_symm_saddle < epsilon_1_saddle) && (epsilon_symm_saddle < epsilon_2_saddle) ) {
728 if (epsilon_1_saddle < epsilon_2_saddle) {
744 eexc_max =
max(eexc1_saddle, eexc2_saddle);
745 eexc_max =
max(eexc_max, e);
747 if ((
G4int)(z) % 2 == 0) {
748 r_e_o_max = 0.3 * (1.0 - 0.2 * (std::pow(z, 2)/a - std::pow(92.0, 2)/238.0));
749 e_pair = 2.0 * 12.0 / std::sqrt(a);
750 if(eexc_max > (e_crit + e_pair)) {
753 if(eexc_max < e_pair) {
756 r_e_o = std::pow((eexc_max - e_crit - e_pair)/e_crit, 2) * r_e_o_max;
782 }
else if (rmode < (ysymm + yasymm1)) {
794 }
else if (imode == 2) {
797 }
else if (imode == 3) {
814 while ( (z1<5.0) || (z2<5.0) ) {
842 re1 =
umass(z1,n1ucd,0.6) +
umass(z2,n2ucd,0.6) +
ecoul(z1,n1ucd,0.6,z2,n2ucd,0.6,2.0);
843 re2 =
umass(z1,n1ucd+1.,0.6) +
umass(z2,n2ucd-1.,0.6) +
ecoul(z1,n1ucd+1.,0.6,z2,n2ucd-1.,0.6,2.0);
844 re3 =
umass(z1,n1ucd+2.,0.6) +
umass(z2,n2ucd-2.,0.6) +
ecoul(z1,n1ucd+2.,0.6,z2,n2ucd-2.,0.6,2.0);
845 reps2 = (re1-2.0*re2+re3) / 2.0;
846 reps1 = re2 - re1 - reps2;
847 rn1_pol = - reps1 / (2.0 * reps2);
848 n1mean = n1ucd + rn1_pol;
850 n1mean = (z1 + 0.5) * n/
z;
860 while (n1r < 5 || n2r < 5) {
864 i_inter = int(n1r + 0.5);
865 n1r = double(i_inter);
874 ne_min = 0.095e0 * a - 20.4e0;
875 if (ne_min < 0) ne_min = 0.0;
876 ne_min = ne_min + e / 8.e0;
885 }
else if (a1 >= 80 && a1 < 110) {
886 ed1_low = (a1-80.)*20./30.;
887 }
else if (a1 >= 110 && a1 < 130) {
888 ed1_low = -(a1-110.)*20./20. + 20.;
889 }
else if (a1 >= 130) {
890 ed1_low = (a1-130.)*20./30.;
895 }
else if (a2 >= 80 && a2 < 110) {
896 ed2_low = (a2-80.)*20./30.;
897 }
else if (a2 >= 110 && a2 < 130) {
898 ed2_low = -(a2-110.)*20./20. + 20.;
899 }
else if (a2 >= 130) {
900 ed2_low = (a2-130.)*20./30.;
903 ed1_high = 20.0*a1/(a1+
a2);
904 ed2_high = 20.0 - ed1_high;
905 ed1 = ed1_low*std::exp(-e/el) + ed1_high*(1-std::exp(-e/el));
906 ed2 = ed2_low*std::exp(-e/el) + ed2_high*(1-std::exp(-e/el));
910 e1 = e*a1/(a1+
a2) + ed1;
911 e2 = e - e*a1/(a1+
a2) + ed2;
922 if ( (icz == -1) || (a1 < 0.0) || (a2 < 0.0) ) {
936 asymm = nsymm + zsymm;
941 cz_symm = 8.0 / std::pow(z,2) * masscurv;
943 wzsymm = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cz_symm) ) ;
957 while ( (z1 < 5.0) || (z2 < 5.0) ) {
961 z1 =
gausshaz(kkk, z1mean, z1width);
977 cn =
umass(zsymm, nsymm+1.0, 0.0) +
umass(zsymm, nsymm-1.0, 0.0)
978 + 1.44 * std::pow(zsymm, 2)/
979 (std::pow(r_null, 2) * std::pow(std::pow(asymm+1.0, 1.0/3.0) + std::pow(asymm-1.0, 1.0/3.0), 2))
980 - 2.0 *
umass(zsymm, nsymm, 0.0) - 1.44 * std::pow(zsymm, 2) /
981 std::pow(r_null * 2.0 *std::pow(asymm, 1.0/3.0), 2);
984 n1width = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cn) + 1.35);
993 re1 =
umass(z1,n1ucd, 0.6) +
umass(z2,n2ucd, 0.6) +
ecoul(z1,n1ucd, 0.6,z2,n2ucd, 0.6,2.0);
994 re2 =
umass(z1,n1ucd+1.,0.6) +
umass(z2,n2ucd-1.,0.6) +
ecoul(z1,n1ucd+1.,0.6,z2,n2ucd-1.,0.6,2.0);
995 re3 =
umass(z1,n1ucd+2.,0.6) +
umass(z2,n2ucd-2.,0.6) +
ecoul(z1,n1ucd+2.,0.6,z2,n2ucd-2.,0.6,2.0);
996 reps2 = (re1-2.0*re2+re3) / 2.0;
997 reps1 = re2 - re1 - reps2;
998 rn1_pol = - reps1 / (2.0 * reps2);
999 n1mean = n1ucd + rn1_pol;
1005 n1r = (float)( (
int)(
gausshaz(k, n1mean,n1width)) );
1012 e2 = e - e*a1/(a1+
a2);
1033 const G4int pSize = 110;
1053 ix = int(y * 100 + 43543000);
1054 if(
mod(ix,2) == 0) {
1066 for(
G4int iRandom = 0; iRandom < pSize; iRandom++) {
1095 G4int count = 0, maxCount = 500;
1097 v1 = 2.0*
haz(k) - 1.0;
1098 v2 = 2.0*
haz(k) - 1.0;
1099 r = std::pow(v1,2) + std::pow(v2,2);
1100 if(++count>maxCount)
1104 fac = std::sqrt(-2.*std::log(r)/r);
1106 fgausshaz = v2*fac*sig+xmoy;
1110 fgausshaz=gset*sig+xmoy;
1162 fractpart = std::modf(number, &intpart);
1167 if(fractpart < 0.5) {
1168 return int(std::floor(number));
1171 return int(std::ceil(number));
1175 if(fractpart < -0.5) {
1176 return int(std::floor(number));
1179 return int(std::ceil(number));
1183 return int(std::floor(number));
1192 mylocaltime = localtime(&mytime);
1195 return(mylocaltime->tm_hour*60*60 + mylocaltime->tm_min*60 + mylocaltime->tm_sec);
1205 return (a - (a/b)*b);
1215 return (a - (a/b)*b);
1227 value = double(std::ceil(a));
1230 value = double(std::floor(a));
1241 value = int(std::ceil(a));
1244 value = int(std::floor(a));
1252 if(a < b && a < c) {
1255 if(b < a && b < c) {
1258 if(c < a && c < b) {
G4int mod(G4int a, G4int b)
G4double dint(G4double a)
G4double umass(G4double z, G4double n, G4double beta)
G4int max(G4int a, G4int b)
double A(double temperature)
G4double dmin1(G4double a, G4double b, G4double c)
void doFission(G4double &A, G4double &Z, G4double &E, G4double &A1, G4double &Z1, G4double &E1, G4double &K1, G4double &A2, G4double &Z2, G4double &E2, G4double &K2)
G4double gausshaz(int k, double xmoy, double sig)
G4double ecoul(G4double z1, G4double n1, G4double beta1, G4double z2, G4double n2, G4double beta2, G4double d)
void even_odd(G4double r_origin, G4double r_even_odd, G4int &i_out)
void fissionDistri(G4double &a, G4double &z, G4double &e, G4double &a1, G4double &z1, G4double &e1, G4double &v1, G4double &a2, G4double &z2, G4double &e2, G4double &v2)
const G4double x[NPOINTSGL]
G4double dmod(G4double a, G4double b)
static const G4double fac
G4int nint(G4double number)
G4int min(G4int a, G4int b)
static const G4double alpha
G4double utilabs(G4double a)