34 #define ABLAXX_IN_GEANT4_MODE 1
44 #ifdef ABLAXX_IN_GEANT4_MODE
50 #ifndef ABLAXX_IN_GEANT4_MODE
108 G4int ZFP1 = 0, AFP1 = 0, AFPIMF = 0, ZFPIMF = 0, ZFP2 = 0, AFP2 = 0;
109 G4double vx_eva = 0.0, vy_eva = 0.0, vz_eva = 0.0;
110 G4double VX_PREF=0.,VY_PREF=0.,VZ_PREF=00,VP1X,VP1Y,VP1Z,VXOUT,VYOUT,VZOUT,V_CM[3],VFP1_CM[3],VFP2_CM[3],VIMF_CM[3],VX2OUT,VY2OUT,VZ2OUT;
111 G4double zf = 0.0, af = 0.0, mtota = 0.0, tkeimf = 0.0, jprf0=0.;
112 G4int ff = 0,afpnew=0,zfpnew=0,aprfp=0,zprfp=0,IOUNSTABLE=0,ILOOP=0,IEV_TAB=0,IEV_TAB_TEMP=0;
113 G4int fimf = 0,INMIN=0,INMAX=0;
115 G4int inum = eventnumber;
140 G4double T_init=0.,T_diff=0.,a_tilda=0.,a_tilda_BU=0., EE_diff=0., EINCL=0., A_FINAL=0., Z_FINAL=0., E_FINAL=0.;
142 G4double A_diff=0.,ASLOPE1,ASLOPE2,A_ACC,ABU_SLOPE, ABU_SUM=0., AMEM=0., ZMEM=0., EMEM=0., JMEM=0., PX_BU_SUM = 0.0, PY_BU_SUM = 0.0, PZ_BU_SUM = 0.0, ETOT_SUM=0., P_BU_SUM=0., ZBU_SUM=0.,Z_Breakup_sum=0.,A_Breakup,Z_Breakup,N_Breakup,G_SYMM,CZ,Sigma_Z,Z_Breakup_Mean,ZTEMP=0.,ATEMP=0.;
144 G4double ETOT_PRF=0.0,PXPRFP=0.,PYPRFP=0.,PZPRFP=0.,PPRFP=0., VX1_BU=0., VY1_BU=0., VZ1_BU=0., VBU2=0., GAMMA_REL=1.0, Eexc_BU_SUM=0., VX_BU_SUM = 0., VY_BU_SUM =0.,VZ_BU_SUM =0., E_tot_BU=0.,EKIN_BU=0.,ZIMFBU=0., AIMFBU=0., ZFFBU=0., AFFBU=0., AFBU=0., ZFBU=0., EEBU=0.,TKEIMFBU=0.,vx_evabu=0.,vy_evabu=0.,vz_evabu=0., Bvalue_BU=0.,P_BU=0.,ETOT_BU=1.,PX_BU=0.,PY_BU=0.,PZ_BU=0.,VX2_BU=0.,VY2_BU=0.,VZ2_BU=0.;
146 G4int ABU_DIFF,ZBU_DIFF,NBU_DIFF;
147 G4int INEWLOOP = 0, ILOOPBU=0;
149 G4double BU_TAB_TEMP[200][5], BU_TAB_TEMP1[200][5];
150 G4double EV_TAB_TEMP[200][5],EV_TEMP[200][5];
151 G4int IMEM_BU[201], IMEM=0;
153 for(
G4int j=0;j<3;j++){
160 for(
G4int I1=0;I1<200;I1++){
161 for(
G4int I2 = 0;I2<11;I2++)
163 for(
G4int I2 = 0;I2<5;I2++){
164 BU_TAB_TEMP[I1][I2] = 0.0;
165 BU_TAB_TEMP1[I1][I2] = 0.0;
166 EV_TAB_TEMP[I1][I2] = 0.0;
169 EV_TEMP[I1][I2] = 0.0;
196 G4double pincl = std::sqrt(pxrem*pxrem + pyrem*pyrem + pzrem*pzrem);
198 G4double ETOT_incl = std::sqrt(pincl*pincl + (AAINCL * amu)*(AAINCL * amu));
199 G4double VX_incl = C * pxrem / ETOT_incl;
200 G4double VY_incl = C * pyrem / ETOT_incl;
201 G4double VZ_incl = C * pzrem / ETOT_incl;
238 a_tilda =
ald->
av*aprf +
ald->
as*std::pow(aprf,2.0/3.0) +
ald->
ak*std::pow(aprf,1.0/3.0);
240 T_init = std::sqrt(EINCL/a_tilda);
244 if(T_diff>0.1 && zprf>2.){
249 for(
G4int i=0;i<5;i++){
256 A_diff =
dint(EE_diff / (8.0 * 5.0 / T_freeze_out));
258 if(A_diff>AAINCL) A_diff = AAINCL;
260 A_FINAL = AAINCL - A_diff;
262 a_tilda =
ald->
av*A_FINAL +
ald->
as*std::pow(A_FINAL,2.0/3.0) +
ald->
ak*std::pow(A_FINAL,1.0/3.0);
266 EE_diff = EINCL - E_FINAL;
278 Z_FINAL =
dint(zprf * A_FINAL / (aprf));
280 if(E_FINAL<0.0) E_FINAL = 0.0;
286 A_diff = AAINCL - aprf;
295 }
else if(A_diff>1.0){
303 a_tilda =
ald->
av*AAINCL +
ald->
as*std::pow(AAINCL,2.0/3.0) +
ald->
ak*std::pow(AAINCL,1.0/3.0);
307 ABU_SLOPE = (ASLOPE1-ASLOPE2)/4.0*(E_FINAL/AAINCL)+
308 ASLOPE1-(ASLOPE1-ASLOPE2)*7.0/4.0;
320 if(ABU_SLOPE > -1.01) ABU_SLOPE = -1.01;
323 Z_Breakup_sum = Z_FINAL;
327 for(
G4int i=0;i<100;i++){
334 std::cout <<
"WARNING: IPOWERLIMHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING A_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED: " << A_Breakup << std::endl;
338 if(A_Breakup>AAINCL)
goto mult4326;
341 std::cout <<
"A_BREAKUP <= 0 " << std::endl;
345 A_ACC = A_ACC + A_Breakup;
349 Z_Breakup_Mean =
dint(A_Breakup * ZAINCL / AAINCL);
351 Z_Breakup_sum = Z_Breakup_sum + Z_Breakup_Mean;
354 G_SYMM = 34.2281 - 5.14037 * E_FINAL/AAINCL;
355 if(E_FINAL/AAINCL < 2.0) G_SYMM = 25.0;
356 if(E_FINAL/AAINCL > 4.0) G_SYMM = 15.0;
361 CZ = 2.0 * G_SYMM * 4.0 / A_Breakup;
365 Sigma_Z = std::sqrt(T_freeze_out/CZ);
373 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING Z_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED: " << A_Breakup <<
" " << Z_Breakup << std::endl;
377 if(Z_Breakup<0.0 )
goto mult4333;
378 if((A_Breakup-Z_Breakup)<0.0)
goto mult4333;
379 if((A_Breakup-Z_Breakup)==0.0 && Z_Breakup!=1.0)
goto mult4333;
381 if(Z_Breakup>=ZAINCL){
384 std::cout <<
"Z_BREAKUP RESAMPLED MORE THAN 10 TIMES; EVENT WILL BE RESAMPLED AGAIN " << std::endl;
394 if(
idnint(A_Breakup-Z_Breakup)<INMIN ||
idnint(A_Breakup-Z_Breakup)>(INMAX+5)){
406 N_Breakup = A_Breakup - Z_Breakup;
409 ABU_SUM = ABU_SUM +
BU_TAB[i][1];
410 ZBU_SUM = ZBU_SUM + BU_TAB[i][0];
413 BU_TAB[I_Breakup][3] = 0.0;
414 I_Breakup = I_Breakup + 1;
415 IMULTBU = IMULTBU + 1;
432 ABU_DIFF =
idnint(ABU_SUM+aprf-AAINCL);
433 ZBU_DIFF =
idnint(ZBU_SUM+zprf-ZAINCL);
434 NBU_DIFF =
idnint((ABU_SUM-ZBU_SUM)+(aprf-zprf)-(AAINCL-ZAINCL));
437 std::cout <<
"WARNING - MORE THAN 200 BU " << IMULTBU << std::endl;
440 std::cout <<
"WARNING - LESS THAN 1 BU " << IMULTBU << std::endl;
444 for(
G4int i=0;i<IMULTBU;i++)
447 while(NBU_DIFF!=0 && ZBU_DIFF!=0){
456 std::cout <<
"WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED." << std::endl;
460 if(IMEM_BU[IEL]==1)
goto mult5555;
461 if(IEL>200)std::cout <<
"5555:" << IEL << RHAZ << IMULTBU << std::endl;
462 if(IEL<0)std::cout <<
"5555:"<< IEL << RHAZ << IMULTBU << std::endl;
465 }
else if(IEL>IMULTBU){
466 N_Breakup =
dint(aprf - zprf -
DSIGN(1.0,
double(NBU_DIFF)));
474 }
else if(IEL>IMULTBU){
475 ZTEMP =
dint(zprf -
DSIGN(1.0,
double(ZBU_DIFF)));
481 if(ZTEMP<1.0 && N_Breakup<1.0){
495 }
else if(IEL>IMULTBU){
497 aprf =
dint(ZTEMP + N_Breakup);
499 NBU_DIFF = NBU_DIFF -
ISIGN(1,NBU_DIFF);
500 ZBU_DIFF = ZBU_DIFF -
ISIGN(1,ZBU_DIFF);
504 for(
G4int i=0;i<IMULTBU;i++)
507 if(NBU_DIFF != 0 && ZBU_DIFF == 0){
508 while(NBU_DIFF > 0 || NBU_DIFF < 0){
514 std::cout <<
"WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED." << std::endl;
518 if(IMEM_BU[IEL]==1)
goto mult5556;
520 if(IPROBA>IMULTBU+1 && NBU_DIFF>0){
521 std::cout <<
"###',IPROBA,IMULTBU,NBU_DIFF,ZBU_DIFF,T_freeze_out" << std::endl;
525 }
else{
if(IEL>IMULTBU)
530 if(IEL>200)std::cout <<
"5556:" << IEL << RHAZ << IMULTBU << std::endl;
531 if(IEL<0)std::cout <<
"5556:"<< IEL << RHAZ << IMULTBU << std::endl;
534 }
else if(IEL>IMULTBU){
543 }
else if(IEL>IMULTBU){
544 ATEMP =
dint(zprf + N_Breakup);
546 if((ATEMP - N_Breakup)<1.0 && N_Breakup<1.0){
558 aprf =
dint(zprf + N_Breakup);
560 NBU_DIFF = NBU_DIFF -
ISIGN(1,NBU_DIFF);
564 for(
G4int i=0;i<IMULTBU;i++)
568 if(ZBU_DIFF != 0 && NBU_DIFF == 0){
569 while(ZBU_DIFF > 0 || ZBU_DIFF < 0){
575 std::cout <<
"WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED." << std::endl;
579 if(IMEM_BU[IEL]==1)
goto mult5557;
581 if(IPROBA>IMULTBU+1 && ZBU_DIFF>0){
582 std::cout <<
"###',IPROBA,IMULTBU,NBU_DIFF,ZBU_DIFF,T_freeze_out" << std::endl;
590 N_Breakup = aprf - zprf;
592 aprf =
dint(zprf + N_Breakup);
597 if(IEL>200)std::cout <<
"5557:" << IEL << RHAZ << IMULTBU << std::endl;
598 if(IEL<0)std::cout <<
"5557:"<< IEL << RHAZ << IMULTBU << std::endl;
602 }
else if(IEL>IMULTBU){
603 N_Breakup =
dint(aprf - zprf);
606 ATEMP =
dint(ZTEMP + N_Breakup);
611 if((ATEMP-ZTEMP)<0.0){
615 if((ATEMP-ZTEMP)<1.0 && ZTEMP<1.0){
625 aprf =
dint(ZTEMP + N_Breakup);
628 ZBU_DIFF = ZBU_DIFF -
ISIGN(1,ZBU_DIFF);
638 for(
G4int i =0;i<IMULTBU;i++){
675 for(
G4int i = 0;i<IMULTBU;i++){
676 ABU_SUM = ABU_SUM +
BU_TAB[i][1];
677 ZBU_SUM = ZBU_SUM + BU_TAB[i][0];
679 ABU_DIFF =
idnint(ABU_SUM-AAINCL);
680 ZBU_DIFF =
idnint(ZBU_SUM-ZAINCL);
682 if(ABU_DIFF!=0 || ZBU_DIFF!=0)
683 std::cout <<
"Problem of mass in BU " << ABU_DIFF <<
" " << ZBU_DIFF << std::endl;
691 AMOMENT(AAINCL,aprf,1,&PXPRFP,&PYPRFP,&PZPRFP);
692 PPRFP = std::sqrt(PXPRFP*PXPRFP + PYPRFP*PYPRFP + PZPRFP*PZPRFP);
695 ETOT_PRF = std::sqrt(PPRFP*PPRFP + (aprf * amu)*(aprf * amu));
696 VX_PREF = C * PXPRFP / ETOT_PRF;
697 VY_PREF = C * PYPRFP / ETOT_PRF;
698 VZ_PREF = C * PZPRFP / ETOT_PRF;
701 tke_bu(zprf,aprf,ZAINCL,AAINCL,&VX1_BU,&VY1_BU,&VZ1_BU);
709 VX_PREF,VY_PREF,VZ_PREF,
710 &VXOUT,&VYOUT,&VZOUT);
717 VBU2 = VX_PREF*VX_PREF + VY_PREF*VY_PREF + VZ_PREF*VZ_PREF;
718 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C*C));
719 ETOT_PRF = aprf * amu / GAMMA_REL;
720 PXPRFP = ETOT_PRF * VX_PREF /
C;
721 PYPRFP = ETOT_PRF * VY_PREF /
C;
722 PZPRFP = ETOT_PRF * VZ_PREF /
C;
736 for(I_Breakup=0;I_Breakup<IMULTBU;I_Breakup++){
739 Eexc_BU_SUM = Eexc_BU_SUM +
BU_TAB[I_Breakup][2];
741 AMOMENT(AAINCL,BU_TAB[I_Breakup][1],1,&PX_BU,&PY_BU,&PZ_BU);
742 P_BU = std::sqrt(PX_BU*PX_BU + PY_BU*PY_BU + PZ_BU*PZ_BU);
745 ETOT_BU = std::sqrt(P_BU*P_BU + (BU_TAB[I_Breakup][1]*amu)*(BU_TAB[I_Breakup][1]*amu));
746 BU_TAB[I_Breakup][4] = C * PX_BU / ETOT_BU;
747 BU_TAB[I_Breakup][5] = C * PY_BU / ETOT_BU;
748 BU_TAB[I_Breakup][6] = C * PZ_BU / ETOT_BU;
750 tke_bu(BU_TAB[I_Breakup][0],BU_TAB[I_Breakup][1],ZAINCL,AAINCL,&VX2_BU,&VY2_BU,&VZ2_BU);
757 BU_TAB[I_Breakup][4],BU_TAB[I_Breakup][5],BU_TAB[I_Breakup][6],
758 &VXOUT,&VYOUT,&VZOUT);
760 BU_TAB[I_Breakup][4] = VXOUT;
761 BU_TAB[I_Breakup][5] = VYOUT;
762 BU_TAB[I_Breakup][6] = VZOUT;
765 VBU2 = BU_TAB[I_Breakup][4]*BU_TAB[I_Breakup][4] +
766 BU_TAB[I_Breakup][5]*BU_TAB[I_Breakup][5] +
767 BU_TAB[I_Breakup][6]*BU_TAB[I_Breakup][6];
768 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C*C));
769 ETOT_BU = BU_TAB[I_Breakup][1]*amu/GAMMA_REL;
770 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
771 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
772 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
774 PX_BU_SUM = PX_BU_SUM + PX_BU;
775 PY_BU_SUM = PY_BU_SUM + PY_BU;
776 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
781 P_BU_SUM = std::sqrt(PX_BU_SUM*PX_BU_SUM + PY_BU_SUM*PY_BU_SUM +
782 PZ_BU_SUM*PZ_BU_SUM);
785 ETOT_SUM = std::sqrt(P_BU_SUM*P_BU_SUM +
786 (AAINCL * amu)*(AAINCL * amu));
788 VX_BU_SUM = C * PX_BU_SUM / ETOT_SUM;
789 VY_BU_SUM = C * PY_BU_SUM / ETOT_SUM;
790 VZ_BU_SUM = C * PZ_BU_SUM / ETOT_SUM;
798 VX_PREF,VY_PREF,VZ_PREF,
799 &VXOUT,&VYOUT,&VZOUT);
805 VBU2 = VX_PREF*VX_PREF + VY_PREF*VY_PREF + VZ_PREF*VZ_PREF;
806 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C*C));
807 ETOT_PRF = aprf * amu / GAMMA_REL;
808 PXPRFP = ETOT_PRF * VX_PREF /
C;
809 PYPRFP = ETOT_PRF * VY_PREF /
C;
810 PZPRFP = ETOT_PRF * VZ_PREF /
C;
821 EKIN_BU = aprf * amu / GAMMA_REL - aprf *
amu;
823 for(I_Breakup=0;I_Breakup<IMULTBU;I_Breakup++){
831 &VXOUT,&VYOUT,&VZOUT);
833 BU_TAB[I_Breakup][4] = VXOUT;
834 BU_TAB[I_Breakup][5] = VYOUT;
835 BU_TAB[I_Breakup][6] = VZOUT;
840 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C*C));
842 ETOT_BU = BU_TAB[I_Breakup][1]*amu/GAMMA_REL;
844 EKIN_BU = EKIN_BU + BU_TAB[I_Breakup][1] * amu /
845 GAMMA_REL - BU_TAB[I_Breakup][1] *
amu;
847 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
848 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
849 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
850 E_tot_BU = E_tot_BU + ETOT_BU;
852 PX_BU_SUM = PX_BU_SUM + PX_BU;
853 PY_BU_SUM = PY_BU_SUM + PY_BU;
854 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
857 if(std::abs(PX_BU_SUM)>10. || std::abs(PY_BU_SUM)>10. ||
858 std::abs(PZ_BU_SUM)>10.){
861 P_BU_SUM = std::sqrt(PX_BU_SUM*PX_BU_SUM + PY_BU_SUM*PY_BU_SUM +
862 PZ_BU_SUM*PZ_BU_SUM);
865 ETOT_SUM = std::sqrt(P_BU_SUM*P_BU_SUM +
866 (AAINCL * amu)*(AAINCL * amu));
868 VX_BU_SUM = C * PX_BU_SUM / ETOT_SUM;
869 VY_BU_SUM = C * PY_BU_SUM / ETOT_SUM;
870 VZ_BU_SUM = C * PZ_BU_SUM / ETOT_SUM;
878 VX_PREF,VY_PREF,VZ_PREF,
879 &VXOUT,&VYOUT,&VZOUT);
885 VBU2 = VX_PREF*VX_PREF + VY_PREF*VY_PREF + VZ_PREF*VZ_PREF;
886 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C*C));
887 ETOT_PRF = aprf * amu / GAMMA_REL;
888 PXPRFP = ETOT_PRF * VX_PREF /
C;
889 PYPRFP = ETOT_PRF * VY_PREF /
C;
890 PZPRFP = ETOT_PRF * VZ_PREF /
C;
901 EKIN_BU = aprf * amu / GAMMA_REL - aprf *
amu;
903 for(I_Breakup=0;I_Breakup<IMULTBU;I_Breakup++){
911 &VXOUT,&VYOUT,&VZOUT);
913 BU_TAB[I_Breakup][4] = VXOUT;
914 BU_TAB[I_Breakup][5] = VYOUT;
915 BU_TAB[I_Breakup][6] = VZOUT;
920 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C*C));
922 ETOT_BU = BU_TAB[I_Breakup][1]*amu/GAMMA_REL;
924 EKIN_BU = EKIN_BU + BU_TAB[I_Breakup][1] * amu /
925 GAMMA_REL - BU_TAB[I_Breakup][1] *
amu;
927 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
928 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
929 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
930 E_tot_BU = E_tot_BU + ETOT_BU;
932 PX_BU_SUM = PX_BU_SUM + PX_BU;
933 PY_BU_SUM = PY_BU_SUM + PY_BU;
934 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
942 for(
G4int i=0;i<IMULTBU;i++){
946 &VP1X,&VP1Y,&VP1Z,BU_TAB_TEMP,&ILOOP);
957 for(
G4int IJ=0;IJ<ILOOP;IJ++){
958 BU_TAB[IMULTBU+INEWLOOP+IJ][0] = BU_TAB_TEMP[IJ][0];
959 BU_TAB[IMULTBU+INEWLOOP+IJ][1] = BU_TAB_TEMP[IJ][1];
960 BU_TAB[IMULTBU+INEWLOOP+IJ][4] = BU_TAB_TEMP[IJ][2];
961 BU_TAB[IMULTBU+INEWLOOP+IJ][5] = BU_TAB_TEMP[IJ][3];
962 BU_TAB[IMULTBU+INEWLOOP+IJ][6] = BU_TAB_TEMP[IJ][4];
963 BU_TAB[IMULTBU+INEWLOOP+IJ][2] = 0.0;
964 BU_TAB[IMULTBU+INEWLOOP+IJ][3] = 0.0;
967 INEWLOOP = INEWLOOP + ILOOP;
974 IMULTBU = IMULTBU + INEWLOOP;
981 for(
G4int i=0;i<IMULTBU;i++){
986 evapora(
BU_TAB[i][0],
BU_TAB[i][1],&EEBU,0.0, &ZFBU, &AFBU, &mtota, &vz_evabu, &vx_evabu,&vy_evabu, &ff, &fimf, &ZIMFBU, &AIMFBU,&TKEIMFBU, &jprfbu, &inttype, &inum,EV_TEMP,&IEV_TAB_TEMP);
991 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
992 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
993 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1000 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1001 &VXOUT,&VYOUT,&VZOUT);
1002 EV_TAB[IJ+IEV_TAB][2] = VXOUT;
1003 EV_TAB[IJ+IEV_TAB][3] = VYOUT;
1004 EV_TAB[IJ+IEV_TAB][4] = VZOUT;
1006 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1016 &VXOUT,&VYOUT,&VZOUT);
1035 G4double EkinR1 = TKEIMFBU * AIMFBU / (AFBU+AIMFBU);
1036 G4double EkinR2 = TKEIMFBU * AFBU / (AFBU+AIMFBU);
1037 G4double V1 = std::sqrt(EkinR1/AFBU) * 1.3887;
1038 G4double V2 = std::sqrt(EkinR2/AIMFBU) * 1.3887;
1040 G4double VPERP1 = std::sqrt(V1*V1 - VZ1_IMF*VZ1_IMF);
1042 G4double VX1_IMF = VPERP1 * std::sin(ALPHA1);
1043 G4double VY1_IMF = VPERP1 * std::cos(ALPHA1);
1048 G4double EEIMFP = EEBU * AFBU /(AFBU + AIMFBU);
1049 G4double EEIMF = EEBU * AIMFBU /(AFBU + AIMFBU);
1052 G4double IINERTTOT = 0.40 * 931.490 * 1.160*1.160 *( std::pow(AIMFBU,5.0/3.0) + std::pow(AFBU,5.0/3.0)) + 931.490 * 1.160*1.160*AIMFBU*AFBU/(AIMFBU+AFBU)*(std::pow(AIMFBU,1./3.) + std::pow(AFBU,1./3.))*(std::pow(AIMFBU,1./3.) + std::pow(AFBU,1./3.));
1054 G4double JPRFHEAVY =
BU_TAB[i][9] * 0.4 * 931.49 * 1.16*1.16 * std::pow(AFBU,5.0/3.0) / IINERTTOT;
1055 G4double JPRFLIGHT =
BU_TAB[i][9] * 0.4 * 931.49 * 1.16*1.16 * std::pow(AIMFBU,5.0/3.0) / IINERTTOT;
1064 &VXOUT,&VYOUT,&VZOUT);
1069 G4double vx1ev_imf=0., vy1ev_imf=0., vz1ev_imf=0., zdummy=0., adummy=0., tkedummy=0.,jprf1=0.;
1072 evapora(ZFBU,AFBU,&EEIMFP,JPRFHEAVY, &ZFFBU, &AFFBU, &mtota, &vz1ev_imf, &vx1ev_imf,&vy1ev_imf, &FFBU1, &FIMFBU1, &zdummy, &adummy,&tkedummy, &jprf1, &inttype, &inum,EV_TEMP,&IEV_TAB_TEMP);
1074 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1075 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1076 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1083 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1084 &VXOUT,&VYOUT,&VZOUT);
1085 EV_TAB[IJ+IEV_TAB][2] = VXOUT;
1086 EV_TAB[IJ+IEV_TAB][3] = VYOUT;
1087 EV_TAB[IJ+IEV_TAB][4] = VZOUT;
1089 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1099 &VXOUT,&VYOUT,&VZOUT);
1109 G4double zffimf, affimf,zdummy1, adummy1, tkedummy1, jprf2, vx2ev_imf, vy2ev_imf, vz2ev_imf;
1111 evapora(ZIMFBU,AIMFBU,&EEIMF,JPRFLIGHT, &zffimf, &affimf, &mtota, &vz2ev_imf, &vx2ev_imf,&vy2ev_imf, &FFBU2, &FIMFBU2, &zdummy1, &adummy1,&tkedummy1, &jprf2, &inttype, &inum,EV_TEMP,&IEV_TAB_TEMP);
1113 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1114 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1115 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1122 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1123 &VXOUT,&VYOUT,&VZOUT);
1126 &VX2OUT,&VY2OUT,&VZ2OUT);
1127 EV_TAB[IJ+IEV_TAB][2] = VX2OUT;
1128 EV_TAB[IJ+IEV_TAB][3] = VY2OUT;
1129 EV_TAB[IJ+IEV_TAB][4] = VZ2OUT;
1131 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1142 &VXOUT,&VYOUT,&VZOUT);
1145 &VX2OUT,&VY2OUT,&VZ2OUT);
1146 BU_TAB[IMULTBU+ILOOPBU][4] = VX2OUT;
1147 BU_TAB[IMULTBU+ILOOPBU][5] = VY2OUT;
1148 BU_TAB[IMULTBU+ILOOPBU][6] = VZ2OUT;
1149 ILOOPBU = ILOOPBU + 1;
1165 IMULTBU = IMULTBU + ILOOPBU;
1173 for(
G4int i=0;i<IMULTBU;i++){
1174 ABU_SUM = ABU_SUM +
BU_TAB[i][8];
1175 ZBU_SUM = ZBU_SUM + BU_TAB[i][7];
1177 BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6],
1178 &VP1X,&VP1Y,&VP1Z,BU_TAB_TEMP1,&ILOOP);
1186 ABU_SUM = ABU_SUM +
G4double(afpnew) - BU_TAB[i][8];
1187 ZBU_SUM = ZBU_SUM +
G4double(zfpnew) - BU_TAB[i][7];
1190 BU_TAB[i][4] = VP1X;
1191 BU_TAB[i][5] = VP1Y;
1192 BU_TAB[i][6] = VP1Z;
1195 for(
G4int IJ=0;IJ<ILOOP;IJ++){
1196 BU_TAB[IMULTBU+INEWLOOP+IJ][7] = BU_TAB_TEMP1[IJ][0];
1197 BU_TAB[IMULTBU+INEWLOOP+IJ][8] = BU_TAB_TEMP1[IJ][1];
1198 BU_TAB[IMULTBU+INEWLOOP+IJ][4] = BU_TAB_TEMP1[IJ][2];
1199 BU_TAB[IMULTBU+INEWLOOP+IJ][5] = BU_TAB_TEMP1[IJ][3];
1200 BU_TAB[IMULTBU+INEWLOOP+IJ][6] = BU_TAB_TEMP1[IJ][4];
1201 BU_TAB[IMULTBU+INEWLOOP+IJ][2] = 0.0;
1202 BU_TAB[IMULTBU+INEWLOOP+IJ][3] = 0.0;
1203 BU_TAB[IMULTBU+INEWLOOP+IJ][0] = BU_TAB[i][0];
1204 BU_TAB[IMULTBU+INEWLOOP+IJ][1] = BU_TAB[i][1];
1205 ABU_SUM = ABU_SUM + BU_TAB[IMULTBU+INEWLOOP+IJ][8];
1206 ZBU_SUM = ZBU_SUM + BU_TAB[IMULTBU+INEWLOOP+IJ][7];
1209 INEWLOOP = INEWLOOP + ILOOP;
1214 IMULTBU = IMULTBU + INEWLOOP;
1218 VX_PREF,VY_PREF,VZ_PREF,
1219 &VXOUT,&VYOUT,&VZOUT);
1224 for(
G4int i=0;i<IMULTBU;i++){
1227 &VXOUT,&VYOUT,&VZOUT);
1232 for(
G4int i=0;i<IEV_TAB;i++){
1235 &VXOUT,&VYOUT,&VZOUT);
1242 if(IMULTBU>200)std::cout <<
"IMULTBU>200 " << IMULTBU << std::endl;
1267 if(zprfp<=2 && zprfp<aprfp){
1285 if(zprfp<=2 && zprfp==aprfp){
1287 VX_PREF, VY_PREF, VZ_PREF,
1288 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
1294 for(
G4int I = 0;I<ILOOP;I++){
1295 for(
G4int IJ = 0; IJ<5; IJ++)
1296 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1298 IEV_TAB = IEV_TAB + ILOOP;
1316 VX_PREF, VY_PREF, VZ_PREF,
1317 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
1323 for(
G4int I = 0;I<ILOOP;I++){
1324 for(
G4int IJ = 0; IJ<5; IJ++)
1325 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1327 IEV_TAB = IEV_TAB + ILOOP;
1342 evapora(zprf,aprf,&ee,jprf, &zf, &af, &mtota, &vz_eva, &vx_eva, &vy_eva, &ff, &fimf, &zimf, &aimf,&tkeimf, &jprf0, &inttype, &inum,EV_TEMP,&IEV_TAB_TEMP);
1344 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1345 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1346 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1353 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1354 &VXOUT,&VYOUT,&VZOUT);
1355 EV_TAB[IJ+IEV_TAB][2] = VXOUT;
1356 EV_TAB[IJ+IEV_TAB][3] = VYOUT;
1357 EV_TAB[IJ+IEV_TAB][4] = VZOUT;
1359 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1365 vx_eva,vy_eva,vz_eva,
1366 &VXOUT,&VYOUT,&VZOUT);
1371 if(ff == 0 && fimf == 0){
1380 VFP1_CM[0] = V_CM[0];
1381 VFP1_CM[1] = V_CM[1];
1382 VFP1_CM[2] = V_CM[2];
1383 for(
G4int j=0;j<3;j++){
1389 if(ff == 1 && fimf == 0) ftype = 1;
1390 if(ff == 0 && fimf == 1) ftype = 2;
1402 G4int IEV_TAB_FIS = 0,imode=0;
1404 G4double vx1_fission=0.,vy1_fission=0.,vz1_fission=0.;
1405 G4double vx2_fission=0.,vy2_fission=0.,vz2_fission=0.;
1406 G4double vx_eva_sc=0.,vy_eva_sc=0.,vz_eva_sc=0.;
1409 &vx1_fission,&vy1_fission,&vz1_fission,
1410 &vx2_fission,&vy2_fission,&vz2_fission,
1411 &ZFP1,&AFP1,&ZFP2,&AFP2,&imode,
1412 &vx_eva_sc,&vy_eva_sc,&vz_eva_sc,EV_TEMP,&IEV_TAB_FIS);
1414 for(
G4int IJ = 0; IJ< IEV_TAB_FIS;IJ++){
1415 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1416 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1423 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1424 &VXOUT,&VYOUT,&VZOUT);
1425 EV_TAB[IJ+IEV_TAB][2] = VXOUT;
1426 EV_TAB[IJ+IEV_TAB][3] = VYOUT;
1427 EV_TAB[IJ+IEV_TAB][4] = VZOUT;
1429 IEV_TAB = IEV_TAB + IEV_TAB_FIS;
1444 V_CM[0],V_CM[1],V_CM[2],
1445 &VXOUT,&VYOUT,&VZOUT);
1448 &VX2OUT,&VY2OUT,&VZ2OUT);
1449 VFP1_CM[0] = VX2OUT;
1450 VFP1_CM[1] = VY2OUT;
1451 VFP1_CM[2] = VZ2OUT;
1458 V_CM[0],V_CM[1],V_CM[2],
1459 &VXOUT,&VYOUT,&VZOUT);
1462 &VX2OUT,&VY2OUT,&VZ2OUT);
1463 VFP2_CM[0] = VX2OUT;
1464 VFP2_CM[1] = VY2OUT;
1465 VFP2_CM[2] = VZ2OUT;
1469 }
else if(ftype == 2){
1477 G4double EkinR1 = tkeimf * aimf / (af+aimf);
1478 G4double EkinR2 = tkeimf * af / (af+aimf);
1480 G4double V2 = std::sqrt(EkinR2/aimf) * 1.3887;
1482 G4double VPERP1 = std::sqrt(V1*V1 - VZ1_IMF*VZ1_IMF);
1484 G4double VX1_IMF = VPERP1 * std::sin(ALPHA1);
1485 G4double VY1_IMF = VPERP1 * std::cos(ALPHA1);
1490 G4double EEIMFP = ee * af /(af + aimf);
1491 G4double EEIMF = ee * aimf /(af + aimf);
1494 G4double IINERTTOT = 0.40 * 931.490 * 1.160*1.160 *( std::pow(aimf,5.0/3.0) + std::pow(af,5.0/3.0)) + 931.490 * 1.160*1.160*aimf*af/(aimf+af)*(std::pow(aimf,1./3.) + std::pow(af,1./3.))*(std::pow(aimf,1./3.) + std::pow(af,1./3.));
1496 G4double JPRFHEAVY = jprf0 * 0.4 * 931.49 * 1.16*1.16 * std::pow(af,5.0/3.0) / IINERTTOT;
1497 G4double JPRFLIGHT = jprf0 * 0.4 * 931.49 * 1.16*1.16 * std::pow(aimf,5.0/3.0) / IINERTTOT;
1498 if(af<2.0) std::cout <<
"RN117-4,AF,ZF,EE,JPRFheavy" << std::endl;
1500 G4double vx1ev_imf=0., vy1ev_imf=0., vz1ev_imf=0., zdummy=0., adummy=0., tkedummy=0.,jprf1=0.;
1502 evapora(zf,af,&EEIMFP,JPRFHEAVY, &zff, &aff, &mtota, &vz1ev_imf, &vx1ev_imf,&vy1ev_imf, &FF11, &FIMF11, &zdummy, &adummy,&tkedummy, &jprf1, &inttype, &inum,EV_TEMP,&IEV_TAB_TEMP);
1504 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1505 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1506 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1513 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1514 &VXOUT,&VYOUT,&VZOUT);
1517 &VX2OUT,&VY2OUT,&VZ2OUT);
1518 EV_TAB[IJ+IEV_TAB][2] = VX2OUT;
1519 EV_TAB[IJ+IEV_TAB][3] = VY2OUT;
1520 EV_TAB[IJ+IEV_TAB][4] = VZ2OUT;
1522 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1531 G4double zffimf, affimf,zdummy1, adummy1, tkedummy1,jprf2,vx2ev_imf,vy2ev_imf,
1534 evapora(zimf,aimf,&EEIMF,JPRFLIGHT, &zffimf, &affimf, &mtota, &vz2ev_imf, &vx2ev_imf,&vy2ev_imf, &FF22, &FIMF22, &zdummy1, &adummy1,&tkedummy1, &jprf2, &inttype, &inum,EV_TEMP,&IEV_TAB_TEMP);
1536 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1537 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1538 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1545 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1546 &VXOUT,&VYOUT,&VZOUT);
1549 &VX2OUT,&VY2OUT,&VZ2OUT);
1550 EV_TAB[IJ+IEV_TAB][2] = VX2OUT;
1551 EV_TAB[IJ+IEV_TAB][3] = VY2OUT;
1552 EV_TAB[IJ+IEV_TAB][4] = VZ2OUT;
1554 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1567 V_CM[0],V_CM[1],V_CM[2],
1568 &VXOUT,&VYOUT,&VZOUT);
1571 &VX2OUT,&VY2OUT,&VZ2OUT);
1572 VIMF_CM[0] = VX2OUT;
1573 VIMF_CM[1] = VY2OUT;
1574 VIMF_CM[2] = VZ2OUT;
1580 V_CM[0],V_CM[1],V_CM[2],
1581 &VXOUT,&VYOUT,&VZOUT);
1584 &VX2OUT,&VY2OUT,&VZ2OUT);
1585 VFP1_CM[0] = VX2OUT;
1586 VFP1_CM[1] = VY2OUT;
1587 VFP1_CM[2] = VZ2OUT;
1589 if(FF11==0 && FIMF11==0){
1598 for(
G4int I=0;I<3;I++)
1602 }
else if(FF11==1 && FIMF11==0){
1615 G4int IEV_TAB_FIS = 0,imode=0;
1617 G4double vx1_fission=0.,vy1_fission=0.,vz1_fission=0.;
1618 G4double vx2_fission=0.,vy2_fission=0.,vz2_fission=0.;
1619 G4double vx_eva_sc=0.,vy_eva_sc=0.,vz_eva_sc=0.;
1622 &vx1_fission,&vy1_fission,&vz1_fission,
1623 &vx2_fission,&vy2_fission,&vz2_fission,
1624 &ZFP1,&AFP1,&ZFP2,&AFP2,&imode,
1625 &vx_eva_sc,&vy_eva_sc,&vz_eva_sc,EV_TEMP,&IEV_TAB_FIS);
1627 for(
G4int IJ = 0; IJ< IEV_TAB_FIS;IJ++){
1628 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1629 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1636 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1637 &VXOUT,&VYOUT,&VZOUT);
1638 EV_TAB[IJ+IEV_TAB][2] = VXOUT;
1639 EV_TAB[IJ+IEV_TAB][3] = VYOUT;
1640 EV_TAB[IJ+IEV_TAB][4] = VZOUT;
1642 IEV_TAB = IEV_TAB + IEV_TAB_FIS;
1655 V_CM[0],V_CM[1],V_CM[2],
1656 &VXOUT,&VYOUT,&VZOUT);
1659 &VX2OUT,&VY2OUT,&VZ2OUT);
1661 VX2OUT,VY2OUT,VZ2OUT,
1662 &VXOUT,&VYOUT,&VZOUT);
1665 &VX2OUT,&VY2OUT,&VZ2OUT);
1666 VFP1_CM[0] = VX2OUT;
1667 VFP1_CM[1] = VY2OUT;
1668 VFP1_CM[2] = VZ2OUT;
1678 V_CM[0],V_CM[1],V_CM[2],
1679 &VXOUT,&VYOUT,&VZOUT);
1682 &VX2OUT,&VY2OUT,&VZ2OUT);
1684 VX2OUT,VY2OUT,VZ2OUT,
1685 &VXOUT,&VYOUT,&VZOUT);
1688 &VX2OUT,&VY2OUT,&VZ2OUT);
1689 VFP2_CM[0] = VX2OUT;
1690 VFP2_CM[1] = VY2OUT;
1691 VFP2_CM[2] = VZ2OUT;
1695 }
else if(FF11==0 && FIMF11==1){
1710 EkinR1 = tkeimf * aimf / (af+aimf);
1711 EkinR2 = tkeimf * af / (af+aimf);
1712 V1 = std::sqrt(EkinR1/af) * 1.3887;
1713 V2 = std::sqrt(EkinR2/aimf) * 1.3887;
1715 VPERP1 = std::sqrt(V1*V1 - VZ1_IMFS*VZ1_IMFS);
1717 G4double VX1_IMFS = VPERP1 * std::sin(ALPHA1);
1718 G4double VY1_IMFS = VPERP1 * std::cos(ALPHA1);
1723 EEIMFP = ee * af /(af + aimf);
1724 EEIMF = ee * aimf /(af + aimf);
1727 IINERTTOT = 0.40 * 931.490 * 1.160*1.160 *( std::pow(aimf,5.0/3.0) + std::pow(af,5.0/3.0)) + 931.490 * 1.160*1.160*aimf*af/(aimf+af)*(std::pow(aimf,1./3.) + std::pow(af,1./3.))*(std::pow(aimf,1./3.) + std::pow(af,1./3.));
1729 JPRFHEAVY = jprf1 * 0.4 * 931.49 * 1.16*1.16 * std::pow(af,5.0/3.0) / IINERTTOT;
1730 JPRFLIGHT = jprf1 * 0.4 * 931.49 * 1.16*1.16 * std::pow(aimf,5.0/3.0) / IINERTTOT;
1732 G4double zffs=0.,affs=0.,vx1ev_imfs=0.,vy1ev_imfs=0.,vz1ev_imfs=0.,jprf3=0.;
1734 evapora(zf,af,&EEIMFP,JPRFHEAVY, &zffs, &affs, &mtota, &vz1ev_imfs, &vx1ev_imfs,&vy1ev_imfs, &FF11, &FIMF11, &zdummy, &adummy,&tkedummy, &jprf3, &inttype, &inum,EV_TEMP,&IEV_TAB_TEMP);
1736 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1737 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1738 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1745 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1746 &VXOUT,&VYOUT,&VZOUT);
1749 &VX2OUT,&VY2OUT,&VZ2OUT);
1750 EV_TAB[IJ+IEV_TAB][2] = VX2OUT;
1751 EV_TAB[IJ+IEV_TAB][3] = VY2OUT;
1752 EV_TAB[IJ+IEV_TAB][4] = VZ2OUT;
1754 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1763 G4double zffimfs=0.,affimfs=0.,vx2ev_imfs=0.,vy2ev_imfs=0.,vz2ev_imfs=0.,jprf4=0.;
1765 evapora(zimf,aimf,&EEIMF,JPRFLIGHT, &zffimfs, &affimfs, &mtota, &vz2ev_imfs, &vx2ev_imfs,&vy2ev_imfs, &FF22, &FIMF22, &zdummy1, &adummy1,&tkedummy1, &jprf4, &inttype, &inum,EV_TEMP,&IEV_TAB_TEMP);
1767 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1768 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1769 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1776 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1777 &VXOUT,&VYOUT,&VZOUT);
1780 &VX2OUT,&VY2OUT,&VZ2OUT);
1781 EV_TAB[IJ+IEV_TAB][2] = VX2OUT;
1782 EV_TAB[IJ+IEV_TAB][3] = VY2OUT;
1783 EV_TAB[IJ+IEV_TAB][4] = VZ2OUT;
1785 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1798 V_CM[0],V_CM[1],V_CM[2],
1799 &VXOUT,&VYOUT,&VZOUT);
1802 &VX2OUT,&VY2OUT,&VZ2OUT);
1804 VX2OUT,VY2OUT,VZ2OUT,
1805 &VXOUT,&VYOUT,&VZOUT);
1808 &VX2OUT,&VY2OUT,&VZ2OUT);
1809 VFP1_CM[0] = VX2OUT;
1810 VFP1_CM[1] = VY2OUT;
1811 VFP1_CM[2] = VZ2OUT;
1819 V_CM[0],V_CM[1],V_CM[2],
1820 &VXOUT,&VYOUT,&VZOUT);
1823 &VX2OUT,&VY2OUT,&VZ2OUT);
1825 VX2OUT,VY2OUT,VZ2OUT,
1826 &VXOUT,&VYOUT,&VZOUT);
1829 &VX2OUT,&VY2OUT,&VZ2OUT);
1830 VFP2_CM[0] = VX2OUT;
1831 VFP2_CM[1] = VY2OUT;
1832 VFP2_CM[2] = VZ2OUT;
1837 if(ftype!=1 && ftype!=21){
1843 VFP1_CM[0],VFP1_CM[1],VFP1_CM[2],
1844 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
1852 for(
G4int I = 0;I<ILOOP;I++){
1853 for(
G4int IJ = 0; IJ<5; IJ++)
1854 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1856 IEV_TAB = IEV_TAB + ILOOP;
1863 VIMF_CM[0],VIMF_CM[1],VIMF_CM[2],
1864 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
1872 for(
G4int I = 0;I<ILOOP;I++){
1873 for(
G4int IJ = 0; IJ<5; IJ++)
1874 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1876 IEV_TAB = IEV_TAB + ILOOP;
1883 VFP2_CM[0],VFP2_CM[1],VFP2_CM[2],
1884 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
1892 for(
G4int I = 0;I<ILOOP;I++){
1893 for(
G4int IJ = 0; IJ<5; IJ++)
1894 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1896 IEV_TAB = IEV_TAB + ILOOP;
1904 if(ftype==1 || ftype==21){
1909 VFP1_CM[0],VFP1_CM[1],VFP1_CM[2],
1910 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
1918 for(
G4int I = 0;I<ILOOP;I++){
1919 for(
G4int IJ = 0; IJ<5; IJ++)
1920 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1922 IEV_TAB = IEV_TAB + ILOOP;
1928 VFP2_CM[0],VFP2_CM[1],VFP2_CM[2],
1929 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
1937 for(
G4int I = 0;I<ILOOP;I++){
1938 for(
G4int IJ = 0; IJ<5; IJ++)
1939 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1941 IEV_TAB = IEV_TAB + ILOOP;
1948 VIMF_CM[0],VIMF_CM[1],VIMF_CM[2],
1949 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
1957 for(
G4int I = 0;I<ILOOP;I++){
1958 for(
G4int IJ = 0; IJ<5; IJ++)
1959 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1961 IEV_TAB = IEV_TAB + ILOOP;
1967 if((ftype == 1 || ftype == 21) && (AFP2<=0 || AFP1<=0 || ZFP2<=0 || ZFP1<=0)){
1968 std::cout <<
"ZFP1:" << ZFP1 << std::endl;
1969 std::cout <<
"AFP1:" << AFP1 << std::endl;
1970 std::cout <<
"ZFP2:" << ZFP2 << std::endl;
1971 std::cout <<
"AFP2:" << AFP2 << std::endl;
1975 EV_TAB[IEV_TAB][0] = ZFP1;
1976 EV_TAB[IEV_TAB][1] = AFP1;
1977 EV_TAB[IEV_TAB][2] = VFP1_CM[0];
1978 EV_TAB[IEV_TAB][3] = VFP1_CM[1];
1979 EV_TAB[IEV_TAB][4] = VFP1_CM[2];
1980 IEV_TAB = IEV_TAB + 1;
1983 EV_TAB[IEV_TAB][0] = ZFP2;
1984 EV_TAB[IEV_TAB][1] = AFP2;
1985 EV_TAB[IEV_TAB][2] = VFP2_CM[0];
1986 EV_TAB[IEV_TAB][3] = VFP2_CM[1];
1987 EV_TAB[IEV_TAB][4] = VFP2_CM[2];
1988 IEV_TAB = IEV_TAB + 1;
1992 EV_TAB[IEV_TAB][0] = ZFPIMF;
1993 EV_TAB[IEV_TAB][1] = AFPIMF;
1994 EV_TAB[IEV_TAB][2] = VIMF_CM[0];
1995 EV_TAB[IEV_TAB][3] = VIMF_CM[1];
1996 EV_TAB[IEV_TAB][4] = VIMF_CM[2];
1997 IEV_TAB = IEV_TAB + 1;
2064 #ifdef ABLAXX_IN_GEANT4_MODE
2069 if(dataInterface->
readData() ==
true) {
2107 for(
G4int i=0;i<13;i++){
2108 for(
G4int j=0;j<154;j++){
2118 mfrldm[j][i] = MP*i+MN*j+
eflmac(i+j,i,1,0);
2123 for(
G4int i=1;i<13;i++){
2124 for(
G4int j=1;j<154;j++){
2144 e0 = 0.285+11.17*std::pow(j+i,-0.464) -0.390-0.00058*(j+i);
2150 e0 = 22.34*std::pow(j+i,-0.464)-0.235;
2157 if((j==i)&&
mod(j,2)==1&&
mod(i,2)==1){
2158 e0 = e0 - 30.0*(1.0/
G4double(j+i));
2172 delete dataInterface;
2259 G4double xv = 0.0, xs = 0.0, xc = 0.0, xa = 0.0;
2261 if ((a <= 0.01) || (z < 0.01)) {
2266 xs = 17.23*std::pow(a,(2.0/3.0));
2269 xc = 0.7*z*(z-1.0)*std::pow((a-1.0),(-1.e0/3.e0));
2276 xa = 23.6*(std::pow((a-2.0*z),2)/
a);
2277 (*el) = xv+xs+xc+xa;
2305 if ( (a1 <= 0) || (z1 <= 0) || ((a1-z1) <= 0) ) {
2314 (*el) =
eflmac(a1,z1,0,refopt4);
2324 (*el) = (*el) + (12.552-0.1436*z1);
2326 if(n1>145&&n1<=152){
2327 (*el) = (*el) + ((152.4-1.77*z1)+(-0.972+0.0113*z1)*n1);
2349 const G4int alpha2Size = 37;
2351 G4double alpha2[alpha2Size] = {0.0, 2.5464e0, 2.4944e0, 2.4410e0, 2.3915e0, 2.3482e0,
2352 2.3014e0, 2.2479e0, 2.1982e0, 2.1432e0, 2.0807e0, 2.0142e0, 1.9419e0,
2353 1.8714e0, 1.8010e0, 1.7272e0, 1.6473e0, 1.5601e0, 1.4526e0, 1.3164e0,
2354 1.1391e0, 0.9662e0, 0.8295e0, 0.7231e0, 0.6360e0, 0.5615e0, 0.4953e0,
2355 0.4354e0, 0.3799e0, 0.3274e0, 0.2779e0, 0.2298e0, 0.1827e0, 0.1373e0,
2356 0.0901e0, 0.0430e0, 0.0000e0};
2361 v = (x - 0.3)/dx + 1.0;
2372 return(alpha2[index] + (alpha2[index+1] - alpha2[index]) / dx * ( x - (0.3e0 + dx*(index-1))));
2387 G4double aa = 0.0,
zz = 0.0, i = 0.0,z2a,C_S,
R,W,G,G1,G2,A_CC;
2397 fissilityResult = std::pow(
zz,2) / aa /50.8830e0 / (1.0e0 - 1.7826e0 * std::pow(i,2));
2402 fissilityResult = std::pow(
zz,2) / aa * std::pow((49.22e0*(1.e0 - 0.3803e0*std::pow(i,2) - 20.489e0*std::pow(i,4))),(-1));
2407 fissilityResult = std::pow(
zz,2) / aa /(48.e0*(1.e0 - 17.22e0*std::pow(i,4)));
2412 C_S = 21.13 * (1.0 - 2.3*i*i);
2413 R = 1.16 * std::pow(aa,1.0/3.0);
2415 G1 = 1.0 - 15.0/8.0*W+21.0/8.0*W*W*W;
2416 G2 = 1.0 + 9.0/2.0*W + 7.0*W*W + 7.0/2.0*W*W*W;
2417 G = 1.0 - 5.0*W*W*(G1 - 3.0/4.0*G2*std::exp(-2.0/W));
2418 A_CC = 3.0/5.0 * 1.44 * G / 1.16;
2419 fissilityResult = z2a * A_CC/(2.0*C_S);
2422 if (fissilityResult > 1.0) {
2423 fissilityResult = 1.0;
2426 if (fissilityResult < 0.0) {
2427 fissilityResult = 0.0;
2430 return fissilityResult;
2433 void G4Abla::evapora(
G4double zprf,
G4double aprf,
G4double *ee_par,
G4double jprf_par,
G4double *zf_par,
G4double *af_par,
G4double *mtota_par,
G4double *vleva_par,
G4double *vxeva_par,
G4double *vyeva_par,
2445 G4int fimf = (*fimf_par);
2447 G4int inttype = (*inttype_par);
2448 G4int inum = (*inum_par);
2540 G4double epsiln = 0.0, probp = 0.0, probd = 0.0, probt = 0.0, probn = 0.0, probhe = 0.0, proba = 0.0, probg = 0.0, probimf=0.0, ptotl = 0.0,
e = 0.0, tcn = 0.0;
2541 G4double sn = 0.0, sbp = 0.0, sbd = 0.0, sbt = 0.0, sbhe = 0.0, sba = 0.0,
x = 0.0, amoins = 0.0, zmoins = 0.0,
sp = 0.0, sd = 0.0, st = 0.0, she = 0.0, sa = 0.0;
2542 G4double ecn = 0.0, ecp = 0.0, ecd = 0.0, ect = 0.0,eche = 0.0,eca = 0.0, ecg = 0.0,
bp = 0.0, bd = 0.0, bt = 0.0, bhe = 0.0, ba = 0.0;
2543 G4double zimf= 0.0,aimf= 0.0,bimf= 0.0,sbimf= 0.0,timf= 0.0;
2544 G4int itest = 0, sortie=0;
2546 G4double ctet1 = 0.0, stet1 = 0.0, phi1 = 0.0;
2550 G4int fgamma = 0, gammadecay = 0;
2552 G4double jprfn=0.0, jprfp=0.0, jprfd=0.0, jprft=0.0, jprfhe=0.0, jprfa=0.0;
2557 const G4double mu2 = 931.494*931.494;
2562 G4int IEV_TAB_TEMP=0;
2564 for(
G4int I1=0;I1<200;I1++)
2565 for(
G4int I2=0;I2<5;I2++)
2566 EV_TEMP[I1][I2] = 0.0;
2576 if(ee<0.|| zf<3.)
goto evapora100;
2577 direct(zf,af,ee,jprf,&probp,&probd,&probt,&probn,&probhe,&proba,&probg,&probimf,&probf,&ptotl,
2578 &sn,&sbp,&sbd,&sbt,&sbhe,&sba,
2579 &ecn,&ecp,&ecd,&ect,&eche,&eca,&ecg,
2580 &
bp,&bd,&bt,&bhe,&ba,&
sp,&sd,&st,&she,&sa,&ef,&ts1,inttype,inum,itest,&sortie,&tcn,
2581 &jprfn, &jprfp, &jprfd, &jprft, &jprfhe, &jprfa, &tsum);
2585 if(ptotl==0.0)
goto evapora100;
2589 if(
e>1e30)std::cout <<
"ERROR AT THE EXIT OF EVAPORA,E>1.D30,AF="<< af <<
" ZF=" << zf << std::endl;
2596 pc = std::sqrt(std::pow((1.0 + (ecn)/9.3956e2),2.) - 1.0) * 9.3956e2;
2602 else if(probp!=0.0){
2606 pc = std::sqrt(std::pow((1.0 + ecp/9.3827e2),2.) - 1.0) * 9.3827e2;
2612 else if(probd!=0.0){
2616 pc = std::sqrt(std::pow((1.0 + ecd/1.875358e3),2) - 1.0) * 1.875358e3;
2622 else if(probt!=0.0){
2626 pc = std::sqrt(std::pow((1.0 + ect/2.80828e3),2) - 1.0) * 2.80828e3;
2632 else if(probhe!=0.0){
2635 epsiln = she + eche;
2636 pc = std::sqrt(std::pow((1.0 + eche/2.80826e3),2) - 1.0) * 2.80826e3;
2642 else{
if(proba!=0.0){
2646 pc = std::sqrt(std::pow((1.0 + eca/3.72834e3),2) - 1.0) * 3.72834e3;
2667 pc = std::sqrt(std::pow((1.0 + eca/3.72834e3),2) - 1.0) * 3.72834e3;
2675 else if (
x < proba+probhe) {
2679 epsiln = she + eche;
2680 pc = std::sqrt(std::pow((1.0 + eche/2.80826e3),2) - 1.0) * 2.80826e3;
2688 else if (
x < proba+probhe+probt) {
2693 pc = std::sqrt(std::pow((1.0 + ect/2.80828e3),2) - 1.0) * 2.80828e3;
2701 else if (
x < proba+probhe+probt+probd) {
2706 pc = std::sqrt(std::pow((1.0 + ecd/1.875358e3),2) - 1.0) * 1.875358e3;
2714 else if (
x < proba+probhe+probt+probd+probp) {
2719 pc = std::sqrt(std::pow((1.0 + ecp/9.3827e2),2) - 1.0) * 9.3827e2;
2727 else if (
x < proba+probhe+probt+probd+probp+probn) {
2732 pc = std::sqrt(std::pow((1.0 + (ecn)/9.3956e2),2.) - 1.0) * 9.3956e2;
2740 else if (
x < proba+probhe+probt+probd+probp+probn+probg) {
2750 if(probp==0.0 && probn==0.0 && probd==0.0 && probt==0.0 && proba==0.0 && probhe==0.0 && probimf==0.0 && probf==0.0)fgamma = 1;
2755 else if (
x < proba+probhe+probt+probd+probp+probn+probg+probimf) {
2762 imf(af,zf,tcn,ee,&zimf,&aimf,&bimf,&sbimf,&timf,jprf);
2764 if(iloop>100)std::cout <<
"Problem in EVAPORA: IMF called > 100 times" << std::endl;
2765 if(zimf>=(zf-2.0))
goto dir1973;
2771 if(zimf==0.0 || aimf==0.0 || sbimf>ee)std::cout <<
"warning: Look in EVAPORA CALL IMF" << std::endl;
2782 tkeimf=
min(2.0*timf,ee-sbimf);
2785 if(tkeimf<=0.0)
goto dir1235;
2786 if(tkeimf>(ee-sbimf) && timf>0.5)
goto dir1235;
2788 tkeimf = tkeimf + bimf;
2792 epsiln = (sbimf-bimf) + tkeimf;
2819 if (ee <= 0.01)ee = 0.01;
2821 if(gammadecay==1 && ee<(epsiln+0.010)){
2822 epsiln = ee - 0.010;
2827 std::cout <<
"***WARNING epsilon<0***" << std::endl;
2835 if (ee <= 0.01)ee = 0.01;
2836 mtota = mtota + malpha;
2839 if(ff==0 && fimf==0){
2841 EV_TEMP[IEV_TAB_TEMP][0] = zmoins;
2842 EV_TEMP[IEV_TAB_TEMP][1] = amoins;
2844 ctet1 = 2.0*rnd - 1.0;
2845 stet1 = std::sqrt(1.0 - std::pow(ctet1,2));
2847 phi1 = rnd*2.0*3.141592654;
2848 G4double xcv = stet1*std::cos(phi1);
2849 G4double ycv = stet1*std::sin(phi1);
2854 G4double ETOT_LP = std::sqrt(pc*pc + amoins*amoins * mu2);
2855 EV_TEMP[IEV_TAB_TEMP][2] = c * pc * xcv / ETOT_LP;
2856 EV_TEMP[IEV_TAB_TEMP][3] = c * pc * ycv / ETOT_LP;
2857 EV_TEMP[IEV_TAB_TEMP][4] = c * pc * zcv / ETOT_LP;
2860 EV_TEMP[IEV_TAB_TEMP][2] = pc * xcv;
2861 EV_TEMP[IEV_TAB_TEMP][3] = pc * ycv;
2862 EV_TEMP[IEV_TAB_TEMP][4] = pc * zcv;
2864 G4double VXOUT=0.,VYOUT=0.,VZOUT=0.;
2866 EV_TEMP[IEV_TAB_TEMP][2],EV_TEMP[IEV_TAB_TEMP][3],
2867 EV_TEMP[IEV_TAB_TEMP][4],
2868 &VXOUT,&VYOUT,&VZOUT);
2869 EV_TEMP[IEV_TAB_TEMP][2] = VXOUT;
2870 EV_TEMP[IEV_TAB_TEMP][3] = VYOUT;
2871 EV_TEMP[IEV_TAB_TEMP][4] = VZOUT;
2874 G4double v2 = std::pow(EV_TEMP[IEV_TAB_TEMP][2],2.) +
2875 std::pow(EV_TEMP[IEV_TAB_TEMP][3],2.) +
2876 std::pow(EV_TEMP[IEV_TAB_TEMP][4],2.);
2877 G4double gamma = 1.0/std::sqrt(1.0 - v2 / (c*c));
2878 G4double etot_lp = amoins*mu * gamma;
2879 pxeva = pxeva - EV_TEMP[IEV_TAB_TEMP][2] * etot_lp / c;
2880 pyeva = pyeva - EV_TEMP[IEV_TAB_TEMP][3] * etot_lp / c;
2881 pleva = pleva - EV_TEMP[IEV_TAB_TEMP][4] * etot_lp / c;
2884 pxeva = pxeva - EV_TEMP[IEV_TAB_TEMP][2];
2885 pyeva = pyeva - EV_TEMP[IEV_TAB_TEMP][3];
2886 pleva = pleva - EV_TEMP[IEV_TAB_TEMP][4];
2888 G4double pteva = std::sqrt(pxeva*pxeva + pyeva*pyeva);
2890 G4double etot = std::sqrt ( pleva*pleva + pteva*pteva + af*af * mu2 );
2891 vxeva = c * pxeva / etot;
2892 vyeva = c * pyeva / etot;
2893 vleva = c * pleva / etot;
2894 IEV_TAB_TEMP = IEV_TAB_TEMP + 1;
2897 if (zf < 3. || (ff == 1) || (fgamma == 1) || (fimf==1)) {
2909 (*tkeimf_par) = tkeimf;
2910 (*mtota_par) = mtota;
2911 (*vleva_par) = vleva;
2912 (*vxeva_par) = vxeva;
2913 (*vyeva_par) = vyeva;
2916 (*inttype_par) = inttype;
2917 (*iev_tab_temp_par)= IEV_TAB_TEMP;
2922 void G4Abla::direct(
G4double zprf,
G4double a,
G4double ee,
G4double jprf,
G4double *probp_par,
G4double *probd_par,
G4double *probt_par,
G4double *probn_par,
G4double *probhe_par,
G4double *proba_par,
G4double *probg_par,
G4double *probimf_par,
G4double *probf_par,
G4double *ptotl_par,
G4double *sn_par,
G4double *sbp_par,
G4double *sbd_par,
G4double *sbt_par,
G4double *sbhe_par,
G4double *sba_par,
G4double *ecn_par,
G4double *ecp_par,
G4double *ecd_par,
G4double *ect_par,
G4double *eche_par,
G4double *eca_par,
G4double *ecg_par,
G4double *bp_par,
G4double *bd_par,
G4double *bt_par,
G4double *bhe_par,
G4double *ba_par,
G4double *sp_par,
G4double *sd_par,
G4double *st_par,
G4double *she_par,
G4double *sa_par,
G4double *ef_par,
G4double *ts1_par,
G4int,
G4int inum,
G4int itest,
G4int *sortie,
G4double *tcn,
G4double *jprfn_par,
G4double *jprfp_par,
G4double *jprfd_par,
G4double *jprft_par,
G4double *jprfhe_par,
G4double *jprfa_par,
G4double *tsum_par)
3081 G4int choice_fisspart = 0;
3184 mglw(a-1.0,zprf,&ma1z);
3185 mglw(a-1.0,zprf-1.0,&ma1z1);
3186 mglw(a-2.0,zprf-1.0,&ma2z1);
3187 mglw(a-3.0,zprf-1.0,&ma3z1);
3188 mglw(a-3.0,zprf-2.0,&ma3z2);
3189 mglw(a-4.0,zprf-2.0,&ma4z2);
3192 if((a-1.)==3.0 && (zprf-1.0)==2.0) ma1z1=-7.7181660;
3193 if((a-1.)==4.0 && (zprf-1.0)==2.0) ma1z1=-28.295992;
3198 sd = ma2z1 - maz - 2.2246;
3199 st = ma3z1 - maz - 8.481977;
3200 she = ma3z2 - maz - 7.7181660;
3201 sa = ma4z2 - maz - 28.295992;
3205 if (zprf <= 1.0e0 || a <= 1.0e0 || (a-zprf) < 0.0) {
3215 if (zprf <= 1.0e0 || a <= 2.0e0 || (a-zprf) < 1.0) {
3225 if (zprf <= 1.0e0 || a <= 3.0e0 || (a-zprf) < 2.0) {
3235 if (a-4.0<=0.0 || zprf<=2.0 || (a-zprf)<2.0) {
3245 if (a-3.0 <= 0.0 || zprf<=2.0 || (a-zprf)<1.0) {
3259 unbound(sn,sp,sd,st,she,sa,bp,bd,bt,bhe,ba,&probf,&probn,&probp,&probd,&probt,&probhe,&proba,&probimf,&probg,&ecn,&ecp,&ecd,&ect,&eche,&eca);
3270 barfit(k,k+j,il,&sbfis,&segs,&selmax);
3276 ef = ef*(4.5114-2.2687*(a-zprf)/zprf);
3278 ef = ef*(3.3931-1.5338*(a-zprf)/zprf);
3282 if((a-zprf)/zprf>1.52)ef=ef*(1.1222-0.10886*(a-zprf)/zprf)-0.1;
3284 if(k>=94&&k<=98&&j<158){
3286 if(
mod(j,2)==0&&
mod(k,2)==0){
3287 if(k>=94){ef = ef-(11.54108*(a-zprf)/zprf-18.074);}
3290 if(
mod(j,2)==1&&
mod(k,2)==1){
3291 if(k>=95){ef = ef-(14.567*(a-zprf)/zprf-23.266);}
3294 if(
mod(j,2)==0&&
mod(k,2)==1){
3295 if(j>=144){ef = ef-(13.662*(a-zprf)/zprf-21.656);}
3298 if(
mod(j,2)==1&&
mod(k,2)==0){
3299 if(j>=144){ef = ef-(13.662*(a-zprf)/zprf-21.656);}
3310 if (ef < 0.0)ef = 0.0;
3341 iinert = 0.4 * 931.49 * 1.16*1.16 * std::pow(a,5.0/3.0)*(1.0 + 0.5*std::sqrt(5./(4.*pi))*defbet);
3342 erot = jprf * jprf * 197.328 * 197.328 /(2. * iinert);
3345 bsbkbc(a,zprf,&bscn,&bkcn,&bccn);
3348 densniv(a,zprf,ee,0.0,&densg,bshell,bscn,bkcn,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprf,0,&qrcn);
3400 lorb(a,a-1.,jprf,ee-sn,&dlout,&sdlout);
3402 if(IDjprf==1) djprf = 0.0;
3403 jprfn = jprf + djprf;
3404 jprfn =
dint(std::abs(jprfn));
3409 iinert = 0.4 * 931.49 * 1.16*1.16 * std::pow(a-1.,5.0/3.0)*(1.0 + 0.5*std::sqrt(5./(4.*pi))*defbet);
3410 erotn = jprfn * jprfn * 197.328 * 197.328 /(2. * iinert);
3411 bsbkbc(a-1.,zprf,&bs,&bk,&bc);
3414 densniv(a-1.0,zprf,ee,sn,&densn,bshell, bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprfn,0,&qr);
3424 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ_NEUT CALLED MORE THAN 100 TIMES" << std::endl;
3433 if(ecn<=0.0)
goto dir1234;
3453 lorb(a,a-1.,jprf,ee-sbp,&dlout,&sdlout);
3455 if(IDjprf==1) djprf = 0.0;
3456 jprfp = jprf + djprf;
3457 jprfp =
dint(std::abs(jprfp));
3462 iinert = 0.4 * 931.49 * 1.16*1.16 * std::pow(a-1.,5.0/3.0)*(1.0 + 0.5*std::sqrt(5./(4.*pi))*defbet);
3463 erotp = jprfp * jprfp * 197.328 * 197.328 /(2. * iinert);
3465 bsbkbc(a-1.,zprf-1.,&bs,&bk,&bc);
3468 densniv(a-1.0,zprf-1.0,ee,sbp,&densp,bshell,bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprfp,0,&qr);
3478 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3487 if(ecp<=0.0)
goto dir1235;
3490 ecp = 2.0 * pt +
bp;
3504 if ((in >= 2) && (iz >= 2)) {
3508 lorb(a,a-2.,jprf,ee-sbd,&dlout,&sdlout);
3510 if(IDjprf==1) djprf = 0.0;
3511 jprfd = jprf + djprf;
3512 jprfd =
dint(std::abs(jprfd));
3517 iinert = 0.4 * 931.49 * 1.16*1.16 * std::pow(a-2.,5.0/3.0)*(1.0 + 0.5*std::sqrt(5./(4.*pi))*defbet);
3518 erotd = jprfd * jprfd * 197.328 * 197.328 /(2. * iinert);
3520 bsbkbc(a-2.,zprf-1.,&bs,&bk,&bc);
3523 densniv(a-2.0,zprf-1.0e0,ee,sbd,&densd,bshell,bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprfd,0,&qr);
3534 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3543 if(ecd<=0.0)
goto dir1236;
3546 ecd = 2.0 * dt + bd;
3560 if ((in >= 3) && (iz >= 2)) {
3564 lorb(a,a-3.,jprf,ee-sbt,&dlout,&sdlout);
3566 if(IDjprf==1) djprf = 0.0;
3567 jprft = jprf + djprf;
3568 jprft =
dint(std::abs(jprft));
3573 iinert = 0.4 * 931.49 * 1.16*1.16 * std::pow(a-3.,5.0/3.0)*(1.0 + 0.5*std::sqrt(5./(4.*pi))*defbet);
3574 erott = jprft * jprft * 197.328 * 197.328 /(2. * iinert);
3576 bsbkbc(a-3.,zprf-1.,&bs,&bk,&bc);
3579 densniv(a-3.0,zprf-1.0,ee,sbt,&denst,bshell,bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprft,0,&qr);
3590 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3599 if(ect<=0.0)
goto dir1237;
3602 ect = 2.0 * tt + bt;
3616 if ((in >= 3) && (iz >= 3)) {
3620 lorb(a,a-4.,jprf,ee-sba,&dlout,&sdlout);
3622 if(IDjprf==1) djprf = 0.0;
3623 jprfa = jprf + djprf;
3624 jprfa =
dint(std::abs(jprfa));
3629 iinert = 0.4 * 931.49 * 1.16*1.16 * std::pow(a-4.,5.0/3.0)*(1.0 + 0.5*std::sqrt(5./(4.*pi))*defbet);
3630 erota = jprfa * jprfa * 197.328 * 197.328 /(2. * iinert);
3632 bsbkbc(a-4.,zprf-2.,&bs,&bk,&bc);
3635 densniv(a-4.0,zprf-2.0,ee,sba,&densa,bshell,bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprfa,0,&qr);
3646 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3655 if(eca<=0.0)
goto dir1238;
3658 eca = 2.0 * at + ba;
3672 if ((in >= 2) && (iz >= 3)) {
3676 lorb(a,a-3.,jprf,ee-sbhe,&dlout,&sdlout);
3678 if(IDjprf==1) djprf = 0.0;
3679 jprfhe = jprf + djprf;
3680 jprfhe =
dint(std::abs(jprfhe));
3685 iinert = 0.4 * 931.49 * 1.16*1.16 * std::pow(a-3.,5.0/3.0)*(1.0 + 0.5*std::sqrt(5./(4.*pi))*defbet);
3686 erothe = jprfhe * jprfhe * 197.328 * 197.328 /(2. * iinert);
3688 bsbkbc(a-3.,zprf-2.,&bs,&bk,&bc);
3691 densniv(a-3.0,zprf-2.0,ee,sbhe,&denshe,bshell,bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprfhe,0,&qr);
3702 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3711 if(eche<=0.0)
goto dir1239;
3714 eche = 2.0 * het + bhe;
3740 gn =
width(a,zprf,1.0,0.0,nt,0.0,sn,ee-erotn)* densn/densg;
3745 gp =
width(a,zprf,1.0,1.0,pt,bp,sbp,ee-erotp)*densp/densg*
pen(a, 1.0, omegap, pt);
3750 gd =
width(a,zprf,2.0,1.0,dt,bd,sbd,ee-erotd)*densd/densg*
pen(a, 2.0, omegad, dt);
3755 gt =
width(a,zprf,3.0,1.0,tt,bt,sbt,ee-erott)*denst/densg*
pen(a, 3.0, omegat, tt);
3760 ghe =
width(a,zprf,3.0,2.0,het,bhe,sbhe,ee-erothe) * denshe/densg*
pen(a, 3.0, omegahe, het);
3765 ga =
width(a,zprf,4.0,2.0,at,ba,sba,ee-erota) * densa/densg*
pen(a, 4.0, omegaa, at);
3773 G4int izcn=0,incn=0,inmin=0,inmax=0,inmi=0,inma=0;
3776 if(fimf_allowed==0 || zprf<=5.0 || a<=7.0){
3794 inmin =
max(inmin,incn-inma);
3795 inmax =
min(inmax,incn-inmi);
3797 inmax =
max(inmax,inmin);
3799 for(
G4int iaimf=izimf+inmin;iaimf<=izimf+inmax;iaimf++){
3801 if(aimf>=a || zimf>=zprf){
3809 sbimf = maimf+mares-mazz+bimf;
3813 iinert= 0.40 * 931.490 * 1.160*1.160 * std::pow(a,5.0/3.0)*(std::pow(aimf,5.0/3.0) + std::pow(a - aimf,5.0/3.0)) + 931.490 * 1.160*1.160 * aimf * (a-aimf) / a *(std::pow(aimf,1.0/3.0) + std::pow(a - aimf,1.0/3.0))*(std::pow(aimf,1.0/3.0) + std::pow(a - aimf,1.0/3.0));
3815 erot = jprf * jprf * 197.328 * 197.328 /(2.0 * iinert);
3818 if(densg==0.0 || ee < (sbimf + erot)){
3824 densniv(a,zprf,ee,sbimf,&densimf,0.0,bsimf,1.0,&timf,0,0,defbetimf,&ecor,jprf,2,&qr);
3826 imfarg = (sbimf+erotcn-erot)/timf;
3827 if(imfarg > 200.0) imfarg = 200.0;
3837 width_imf =
width(a,zprf,aimf,zimf,timf,bimf,sbimf,ee-erot)*std::exp(-imfarg)*qr/qrcn;
3840 gimf3 = gimf3 + width_imf;
3854 inmin =
max(inmin,incn-inma);
3855 inmax =
min(inmax,incn-inmi);
3857 inmax =
max(inmax,inmin);
3859 for(
G4int iaimf=izimf+inmin;iaimf<=izimf+inmax;iaimf++){
3861 if(aimf>=a || zimf>=zprf){
3869 sbimf = maimf+mares-mazz+bimf;
3873 iinert= 0.40 * 931.490 * 1.160*1.160 * std::pow(a,5.0/3.0)*(std::pow(aimf,5.0/3.0) + std::pow(a - aimf,5.0/3.0)) + 931.490 * 1.160*1.160 * aimf * (a-aimf) / a *(std::pow(aimf,1.0/3.0) + std::pow(a - aimf,1.0/3.0))*(std::pow(aimf,1.0/3.0) + std::pow(a - aimf,1.0/3.0));
3875 erot = jprf * jprf * 197.328 * 197.328 /(2.0 * iinert);
3878 if(densg==0.0 || ee < (sbimf + erot)){
3884 densniv(a,zprf,ee,sbimf,&densimf,0.0,bsimf,1.0,&timf,0,0,defbetimf,&ecor,jprf,2,&qr);
3886 imfarg = (sbimf+erotcn-erot)/timf;
3887 if(imfarg > 200.0) imfarg = 200.0;
3896 width_imf =
width(a,zprf,aimf,zimf,timf,bimf,sbimf,ee-erot)*std::exp(-imfarg)*qr/qrcn;
3899 gimf5 = gimf5 + width_imf;
3904 if(gimf3<=0.0 || gimf5<=0.0){
3910 b_imf = (std::log10(gimf3) - std::log10(gimf5))/(std::log10(3.0)-std::log10(5.0));
3912 if(b_imf >= -1.01) b_imf = -1.01;
3913 if(b_imf <= -100.0) {
3920 a_imf = gimf3 / std::pow(3.0,b_imf);
3921 gimf = a_imf * ( std::pow(zprf,b_imf+1.0) - std::pow(3.0,b_imf+1.0)) /(b_imf + 1.0);
3925 if(gimf < 1.
e-10) gimf = 0.0;
3932 pa = (
ald->
av)*a + (
ald->
as)*std::pow(a,2./3.) + (
ald->
ak)*std::pow(a,1./3.);
3933 gamma = 2.5 * pa * std::pow(a,-4./3.);
3939 gtemp = 17.60/(std::pow(a,0.699) * std::sqrt(gfactor));
3942 gg = 0.624e-9*std::pow(a,1.6)*std::pow(gtemp,5.);
3952 gsum = ga + ghe + gd + gt + gp + gn + gimf + gg;
3962 if(
fiss->
ifis==0 || (zprf*zprf/a<=22.74 && zprf<60.)){
3970 fission_width(zprf,a,ee,bssp,bksp,ef,y,&gf,&temp,jprf,0,1,
fiss->
optcol,
fiss->
optshp,densg);
3989 gtotal = ga + ghe + gp + gd + gt + gn + gg +gimf + gf;
4007 probhe = ghe/gtotal;
4010 probimf = gimf/gtotal;
4023 fomega_sp(a,y,&mfcd,&omegasp,&homegasp);
4027 fomega_gs(a,zprf,&k1,&omegags,&homegags);
4040 part_fiss(
fiss->
bet,gsum,gf,y,tauc,ts1,tsum, &choice_fisspart,zprf,a,ft,&t_lapse,&gff);
4044 tsum = tsum + t_lapse;
4047 if(choice_fisspart==2){
4062 gtotal=ga + ghe + gp + gd + gt + gn + gimf + gg;
4080 probhe = ghe/gtotal;
4083 probimf = gimf/gtotal;
4093 gtotal = ga + ghe + gp + gd + gt + gn + gg + gimf + gf;
4099 probhe = ghe/gtotal;
4102 probimf = gimf/gtotal;
4107 gtotal = ga + ghe + gp + gd + gt + gn + gg + gimf;
4124 probhe = ghe/gtotal;
4127 probimf = gimf/gtotal;
4131 ptotl = probp+probd+probt+probn+probhe+proba+probg+probimf+probf;
4137 (*probp_par) = probp;
4138 (*probd_par) = probd;
4139 (*probt_par) = probt;
4140 (*probn_par) = probn;
4141 (*probhe_par) = probhe;
4142 (*proba_par) = proba;
4143 (*probg_par) = probg;
4144 (*probimf_par) = probimf;
4145 (*probf_par) = probf;
4146 (*ptotl_par) = ptotl;
4172 (*jprfn_par) = jprfn;
4173 (*jprfp_par) = jprfp;
4174 (*jprfd_par) = jprfd;
4175 (*jprft_par) = jprft;
4176 (*jprfhe_par) = jprfhe;
4177 (*jprfa_par) = jprfa;
4182 void G4Abla::densniv(
G4double a,
G4double z,
G4double ee,
G4double esous,
G4double *dens,
G4double bshell,
G4double bsin,
G4double bkin,
G4double *temp,
G4int optshp,
G4int optcol,
G4double defbet,
G4double *ecor,
G4double jprf,
G4int ifis,
G4double *qr)
4277 G4double pi6 = std::pow(3.1415926535,2) / 6.0;
4289 if(afp<=20) BSHELLCT = 0.0;
4322 pa = (
ald->
av)*a + (
ald->
as)*std::pow(a,(2.e0/3.e0)) + (
ald->
ak)*std::pow(a,(1.e0/3.e0));
4324 pa = (
ald->
av)*a + (
ald->
as)*bsin*std::pow(a,(2.e0/3.e0)) + (
ald->
ak)*bkin*std::pow(a,(1.e0/3.e0));
4326 gamma = 2.5 * pa * std::pow(a,-4.0/3.0);
4331 if(ifis==0&&bs!=1.0){
4334 if(ponq>700.0) ponq = 700.0;
4335 bs = 1.0/(1.0+std::exp(-ponq)) + 1.0/(1.0+std::exp(ponq)) * bsin;
4336 bk = 1.0/(1.0+std::exp(-ponq)) + 1.0/(1.0+std::exp(ponq)) * bkin;
4341 pa = (
ald->
av)*a + (
ald->
as)*std::pow(a,(2.e0/3.e0)) + (
ald->
ak)*std::pow(a,(1.e0/3.e0));
4344 pa = (
ald->
av)*a + (
ald->
as)*bs*std::pow(a,(2.e0/3.e0)) + (
ald->
ak)*bk*std::pow(a,(1.e0/3.e0));
4347 gamma = 2.5 * pa * std::pow(a,-4.0/3.0);
4353 ecr = pa*17.60/(std::pow(a,0.699) * std::sqrt(1.0+gamma*BSHELLCT))*17.60/(std::pow(a,0.699) * std::sqrt(1.0+gamma*BSHELLCT));
4373 deltpp = -0.25e0* std::pow((delta0/std::sqrt(a)),2) * pa /pi6 + 22.34e0*std::pow(a,-0.464)-0.235;
4377 e=e-(0.285+11.17*std::pow(a,-0.464)-0.390-0.00058*
a);
4381 e=e-(22.34*std::pow(a,-0.464)-0.235);
4405 ponfe = -2.5*pa*e*std::pow(a,(-4.0/3.0));
4407 if (ponfe < -700.0) {
4410 fe = 1.0 - std::exp(ponfe);
4413 he = 1.0 - std::pow((1.0 - e/ecr),2);
4420 fecor = e + deltau*fe + deltpp*he;
4427 y1 = std::sqrt(pa*fecor);
4428 for(
G4int j = 0; j < 5; j++) {
4429 y2 = pa*fecor*(1.e0-std::exp(-y1));
4434 fdens = std::exp(y0*fecor)/ (std::pow((std::pow(fecor,3)*y0),0.5)*std::pow((1.0-0.5*y0*fecor*std::exp(-y1)),0.5))* std::exp(y1)*(1.0-std::exp(-y1))*0.1477045;
4437 y11 = std::sqrt(pa*ecor1);
4438 for(
G4int j = 0; j < 7; j++) {
4439 y21 = pa*ecor1*(1.0-std::exp(-y11));
4440 y11 = std::sqrt(y21);
4444 fdens = fdens*std::pow((y01/y0),1.5);
4445 ftemp = ftemp*std::pow((y01/y0),1.5);
4449 ponniv = 2.0*std::sqrt(pa*fecor);
4450 if (ponniv > 700.0) {
4454 fdens = 0.1477045 * std::exp(ponniv)/(std::pow(pa,0.25)*std::pow(fecor,1.25));
4455 ftemp = std::sqrt(fecor/pa);
4461 if(IOPTCT==0)
goto densniv100;
4462 tempct = 17.60/( std::pow(a,0.699) * std::sqrt(1.+gamma*BSHELLCT));
4473 if (IPARITE == 1) { e0 = 0.285+11.17*std::pow(a,-0.464) - 0.390-0.00058*
a;}
4475 if (IPARITE == 2) { e0 = 22.34*std::pow(a,-0.464)-0.235;}
4477 if (IPARITE == 0){ e0 = 0.0;}
4479 ponniv = (ein-e0)/tempct;
4480 if(ifis!=1) ponniv =
max(0.0,(ein-e0)/tempct);
4481 if(ponniv>700.0){ ponniv = 700.0;}
4482 densct = std::exp(ponniv)/tempct*std::exp(0.079*BSHELLCT/tempct);
4486 if(elim>=ecr&&densfm<=densct){
4494 if(elim>=ecr&&tfm>=tempct){
4502 ponniv = (ein)/tempct;
4503 if(ponniv>700.0){ ponniv = 700.0;}
4504 densct = std::exp(ponniv)/tempct;
4506 if(ein>=ecr && densfm<=densct){
4516 if(ein>=ecr && tfm>=tempct){
4531 ftemp = 17.60/( std::pow(a,0.699) * std::sqrt(1.0+gamma*BSHELLCT));
4544 fnorm = std::pow(1.16,2)*931.49*1.e-2/(9.0* std::pow(6.582122,2));
4546 if(ifis==0 || ifis==2){
4552 fp_per = 0.4*std::pow(a,5.0/3.0)*fnorm*(1.0+0.50*defbet*std::sqrt(5.0/(4.0*pi)));
4553 fp_par = 0.40*std::pow(a,5.0/3.0)*fnorm*(1.0-defbet*std::sqrt(5.0/(4.0*pi)));
4562 fp_per = 2.0/5.0*std::pow(a,5.0/3.0)*fnorm*(1.0+7.0/6.0*defbet*(1.0+1396.0/255.0*defbet));
4564 fp_par = 2.0/5.0*std::pow(a,5.0/3.0)*fnorm*(1.0-7.0/3.0*defbet*(1.0-389.0/255.0*defbet));
4571 fp_per = 0.4*std::pow(a,5.0/3.0)*fnorm*3.50*(1.0 + std::pow(defbet,5.))/std::pow(1.0 + defbet*defbet*defbet,5.0/3.0);
4572 fp_par = 0.4*std::pow(a,5.0/3.0)*fnorm*(1.0 + std::pow(defbet,5.0))/std::pow(1.0 + defbet*defbet*defbet,5.0/3.0);
4576 if(fp_par<0.0)fp_par=0.0;
4577 if(fp_per<0.0)fp_per=0.0;
4579 sig_per = std::sqrt(fp_per * ftemp);
4580 sig_par = std::sqrt(fp_par * ftemp);
4582 sigma2 = sig_per*sig_per + sig_par*sig_par;
4583 jfact = (2.*jprf+1.)*std::exp(-1.*jprf*(jprf+1.0)/(2.0*sigma2))/(std::sqrt(8.0*3.1415)*std::pow(sigma2,1.5));
4584 erot = jprf*jprf/(2.0*std::sqrt(fp_par*fp_par+fp_per*fp_per));
4588 qrot(z,a,defbet,sig_per,fecor-erot,&fqr);
4594 fdens = fdens * fqr *jfact;
4596 if(fdens<1e-300)fdens=0.0;
4632 G4double ponq = 0.0, dn = 0.0,
n = 0.0, dz = 0.0;
4633 G4int distn,distz,ndist, zdist;
4634 G4int nmn[8]= {2, 8, 14, 20, 28, 50, 82, 126};
4635 G4int nmz[8]= {2, 8, 14, 20, 28, 50, 82, 126};
4639 if(std::abs(bet)<=0.15){
4650 for(
G4int i =0;i<8;i++){
4651 ndist = std::fabs(
idnint(
n) - nmn[i]);
4652 if(ndist < distn) distn = ndist;
4653 zdist = std::fabs(
idnint(z) - nmz[i]);
4654 if(zdist < distz) distz = zdist;
4660 bet = 0.022 + 0.003*dn + 0.002*dz;
4662 sig = 75.0*std::pow(bet,2.) * sig;
4666 ponq = (u - ucr)/dcr;
4674 (*qr) = 1.0/(1.0 + std::exp(ponq)) * (sig - 1.0) + 1.0;
4695 for(
G4int i = 2; i <
n; i++) {
4722 G4double z = 0.0,
n = 0.0,
a = 0.0, av = 0.0, as = 0.0;
4723 G4double a0 = 0.0,
c1 = 0.0, c4 = 0.0, b1 = 0.0, b3 = 0.0;
4725 G4double r0 = 0.0, kf = 0.0, ks = 0.0;
4726 G4double kv = 0.0, rp = 0.0, ay = 0.0, aden = 0.0, x0 = 0.0, y0 = 0.0;
4727 G4double esq = 0.0, ael = 0.0, i = 0.0, e0 = 0.0;
4777 if(flag==1){
goto eflmac311;}
4787 c1 = 3.0/5.0*esq/r0;
4788 c4 = 5.0/4.0*std::pow((3.0/(2.0*pi)),(2.0/3.0)) *
c1;
4789 kf = std::pow((9.0*pi*z/(4.0*
a)),(1.0/3.0))/r0;
4791 ff = -1.0/8.0*rp*rp*esq/std::pow(r0,3) * (145.0/48.0 - 327.0/2880.0*std::pow(kf,2) * std::pow(rp,2) + 1527.0/1209600.0*std::pow(kf,4) * std::pow(rp,4));
4794 x0 = r0 * std::pow(
a,(1.0/3.0)) / ay;
4795 y0 = r0 * std::pow(
a,(1.0/3.0)) / aden;
4797 b1 = 1.0 - 3.0/(std::pow(x0,2)) + (1.0 + x0) * (2.0 + 3.0/x0 + 3.0/std::pow(x0,2)) * std::exp(-2.0*x0);
4799 b3 = 1.0 - 5.0/std::pow(y0,2) * (1.0 - 15.0/(8.0*y0) + 21.0/(8.0 * std::pow(y0,3))
4800 - 3.0/4.0 * (1.0 + 9.0/(2.0*y0) + 7.0/std::pow(y0,2)
4801 + 7.0/(2.0 * std::pow(y0,3))) * std::exp(-2.0*y0));
4805 efl = -1.0 * av*(1.0 - kv*i*i)*
a + as*(1.0 - ks*i*i)*b1 * std::pow(
a,(2.0/3.0)) + a0
4806 +
c1*z*z*b3/std::pow(
a,(1.0/3.0)) - c4*std::pow(z,(4.0/3.0))/std::pow(
a,(1.e0/3.e0))
4807 + ff*std::pow(z,2)/
a -ca*(
n-
z) - ael * std::pow(z,(2.39e0));
4809 efl = efl +
w*std::abs(i);
4814 if (in==iz && (
mod(in,2) == 1) && (
mod(iz,2) == 1) && in>0.) {
4830 e0 = 0.285+11.17*std::pow(
a,-0.464) -0.390-0.00058*(
a);
4836 e0 = 22.34*std::pow(
a,-0.464)-0.235;
4848 return eflmacResult;
4873 (*del) = -12.0/std::sqrt(a);
4876 (*del) = 12.0/std::sqrt(a);
4888 G4double n1 = 0.0, n2 = 0.0, n3 = 0.0;
4924 if (bet/(std::sqrt(2.0)*10.0*(homega/6.582122)) <= 1.0) {
4925 tauResult = std::log(10.0*ef/t)/(bet*1.0e21);
4928 tauResult = std::log(10.0*ef/t)/ (2.0*std::pow((10.0*homega/6.582122),2))*(bet*1.0
e-21);
4940 G4double rel = bet/(20.0*homega/6.582122);
4941 G4double cramResult = std::sqrt(1.0 + std::pow(rel,2)) - rel;
4944 if (cramResult > 1.0) {
4966 const G4int bsbkSize = 54;
4968 G4double bk[bsbkSize] = {0.0, 1.00000,1.00087,1.00352,1.00799,1.01433,1.02265,1.03306,
4969 1.04576,1.06099,1.07910,1.10056,1.12603,1.15651,1.19348,
4970 1.23915,1.29590,1.35951,1.41013,1.44103,1.46026,1.47339,
4971 1.48308,1.49068,1.49692,1.50226,1.50694,1.51114,1.51502,
4972 1.51864,1.52208,1.52539,1.52861,1.53177,1.53490,1.53803,
4973 1.54117,1.54473,1.54762,1.55096,1.55440,1.55798,1.56173,
4974 1.56567,1.56980,1.57413,1.57860,1.58301,1.58688,1.58688,
4975 1.58688,1.58740,1.58740, 0.0};
4977 G4double bs[bsbkSize] = {0.0, 1.00000,1.00086,1.00338,1.00750,1.01319,
4978 1.02044,1.02927,1.03974,
4979 1.05195,1.06604,1.08224,1.10085,1.12229,1.14717,1.17623,1.20963,
4980 1.24296,1.26532,1.27619,1.28126,1.28362,1.28458,1.28477,1.28450,
4981 1.28394,1.28320,1.28235,1.28141,1.28042,1.27941,1.27837,1.27732,
4982 1.27627,1.27522,1.27418,1.27314,1.27210,1.27108,1.27006,1.26906,
4983 1.26806,1.26707,1.26610,1.26514,1.26418,1.26325,1.26233,1.26147,
4984 1.26147,1.26147,1.25992,1.25992, 0.0};
4986 i =
idint(y/(2.0
e-02)) + 1;
4988 if((i + 1) >= bsbkSize) {
4996 bipolResult = bs[i] + (bs[i+1] - bs[i])/2.0
e-02 * (y - 2.0
e-02*(i - 1));
4999 bipolResult = bk[i] + (bk[i+1] - bk[i])/2.0
e-02 * (y - 2.0
e-02*(i - 1));
5014 ES0 = 20.760*std::pow(AF,2.0/3.0);
5017 MR02 = std::pow(AF,5.0/3.0)*1.0340*0.010*1.175*1.175;
5019 (*MFCD) = MR02 * 3.0/10.0*(1.0+3.0*
Y);
5021 OMEGA = std::sqrt(ES0/MR02)*std::sqrt(8.0/3.0*Y*(1.0+304.0*Y/255.0));
5023 HOMEGA = 6.58122*OMEGA/10.0;
5038 G4double OMEGA,HOMEGA,MR02,MINERT,
C,fk1;
5040 MR02 = std::pow(AF,5.0/3.0)*1.0340*0.01*1.175*1.175;
5041 MINERT = 3.*MR02/10.0;
5042 C = 17.9439*(1.-1.7826*std::pow((AF-2.0*ZF)/AF,2));
5043 fk1 = 0.4*C*std::pow(AF,2.0/3.0)-0.1464*std::pow(ZF,2)/std::pow(AF,1./3.);
5044 OMEGA = std::sqrt(fk1/MINERT);
5045 HOMEGA = 6.58122*OMEGA/10.0;
5080 BARR = 1.345 * Z1 * Z2 / RMAX;
5082 OMEGA = 4.5 / 197.3287;
5159 for(
G4int init_i = 0; init_i < 7; init_i++) {
5163 for(
G4int init_i = 0; init_i < 10; init_i++) {
5167 G4double a = 0.0,
z = 0.0, amin = 0.0, amax = 0.0, amin2 = 0.0;
5168 G4double amax2 = 0.0, aa = 0.0,
zz = 0.0, bfis = 0.0;
5169 G4double bfis0 = 0.0, ell = 0.0, el = 0.0, egs = 0.0, el80 = 0.0, el20 = 0.0;
5170 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;
5171 G4double aj = 0.0, ak = 0.0, a1 = 0.0, a2 = 0.0;
5173 G4int i = 0, j = 0, k = 0,
m = 0;
5176 G4double emncof[4][5] = {{-9.01100e+2,-1.40818e+3, 2.77000e+3,-7.06695e+2, 8.89867e+2},
5177 {1.35355e+4,-2.03847e+4, 1.09384e+4,-4.86297e+3,-6.18603e+2},
5178 {-3.26367e+3, 1.62447e+3, 1.36856e+3, 1.31731e+3, 1.53372e+2},
5179 {7.48863e+3,-1.21581e+4, 5.50281e+3,-1.33630e+3, 5.05367e-2}};
5181 G4double elmcof[4][5] = {{1.84542e+3,-5.64002e+3, 5.66730e+3,-3.15150e+3, 9.54160e+2},
5182 {-2.24577e+3, 8.56133e+3,-9.67348e+3, 5.81744e+3,-1.86997e+3},
5183 {2.79772e+3,-8.73073e+3, 9.19706e+3,-4.91900e+3, 1.37283e+3},
5184 {-3.01866e+1, 1.41161e+3,-2.85919e+3, 2.13016e+3,-6.49072e+2}};
5186 G4double emxcof[4][6] = {{9.43596e4,-2.241997e5,2.223237e5,-1.324408e5,4.68922e4,-8.83568e3},
5187 {-1.655827e5,4.062365e5,-4.236128e5,2.66837e5,-9.93242e4,1.90644e4},
5188 {1.705447e5,-4.032e5,3.970312e5,-2.313704e5,7.81147e4,-1.322775e4},
5189 {-9.274555e4,2.278093e5,-2.422225e5,1.55431e5,-5.78742e4,9.97505e3}};
5191 G4double elzcof[7][7] = {{5.11819909e+5,-1.30303186e+6, 1.90119870e+6,-1.20628242e+6, 5.68208488e+5, 5.48346483e+4,-2.45883052e+4},
5192 {-1.13269453e+6, 2.97764590e+6,-4.54326326e+6, 3.00464870e+6, -1.44989274e+6,-1.02026610e+5, 6.27959815e+4},
5193 {1.37543304e+6,-3.65808988e+6, 5.47798999e+6,-3.78109283e+6, 1.84131765e+6, 1.53669695e+4,-6.96817834e+4},
5194 {-8.56559835e+5, 2.48872266e+6,-4.07349128e+6, 3.12835899e+6, -1.62394090e+6, 1.19797378e+5, 4.25737058e+4},
5195 {3.28723311e+5,-1.09892175e+6, 2.03997269e+6,-1.77185718e+6, 9.96051545e+5,-1.53305699e+5,-1.12982954e+4},
5196 {4.15850238e+4, 7.29653408e+4,-4.93776346e+5, 6.01254680e+5, -4.01308292e+5, 9.65968391e+4,-3.49596027e+3},
5197 {-1.82751044e+5, 3.91386300e+5,-3.03639248e+5, 1.15782417e+5, -4.24399280e+3,-6.11477247e+3, 3.66982647e+2}};
5199 const G4int sizex = 5;
5200 const G4int sizey = 6;
5201 const G4int sizez = 4;
5203 G4double egscof[sizey][sizey][sizez];
5205 G4double egs1[sizey][sizex] = {{1.927813e5, 7.666859e5, 6.628436e5, 1.586504e5,-7.786476e3},
5206 {-4.499687e5,-1.784644e6,-1.546968e6,-4.020658e5,-3.929522e3},
5207 {4.667741e5, 1.849838e6, 1.641313e6, 5.229787e5, 5.928137e4},
5208 {-3.017927e5,-1.206483e6,-1.124685e6,-4.478641e5,-8.682323e4},
5209 {1.226517e5, 5.015667e5, 5.032605e5, 2.404477e5, 5.603301e4},
5210 {-1.752824e4,-7.411621e4,-7.989019e4,-4.175486e4,-1.024194e4}};
5212 G4double egs2[sizey][sizex] = {{-6.459162e5,-2.903581e6,-3.048551e6,-1.004411e6,-6.558220e4},
5213 {1.469853e6, 6.564615e6, 6.843078e6, 2.280839e6, 1.802023e5},
5214 {-1.435116e6,-6.322470e6,-6.531834e6,-2.298744e6,-2.639612e5},
5215 {8.665296e5, 3.769159e6, 3.899685e6, 1.520520e6, 2.498728e5},
5216 {-3.302885e5,-1.429313e6,-1.512075e6,-6.744828e5,-1.398771e5},
5217 {4.958167e4, 2.178202e5, 2.400617e5, 1.167815e5, 2.663901e4}};
5219 G4double egs3[sizey][sizex] = {{3.117030e5, 1.195474e6, 9.036289e5, 6.876190e4,-6.814556e4},
5220 {-7.394913e5,-2.826468e6,-2.152757e6,-2.459553e5, 1.101414e5},
5221 {7.918994e5, 3.030439e6, 2.412611e6, 5.228065e5, 8.542465e3},
5222 {-5.421004e5,-2.102672e6,-1.813959e6,-6.251700e5,-1.184348e5},
5223 {2.370771e5, 9.459043e5, 9.026235e5, 4.116799e5, 1.001348e5},
5224 {-4.227664e4,-1.738756e5,-1.795906e5,-9.292141e4,-2.397528e4}};
5226 G4double egs4[sizey][sizex] = {{-1.072763e5,-5.973532e5,-6.151814e5, 7.371898e4, 1.255490e5},
5227 {2.298769e5, 1.265001e6, 1.252798e6,-2.306276e5,-2.845824e5},
5228 {-2.093664e5,-1.100874e6,-1.009313e6, 2.705945e5, 2.506562e5},
5229 {1.274613e5, 6.190307e5, 5.262822e5,-1.336039e5,-1.115865e5},
5230 {-5.715764e4,-2.560989e5,-2.228781e5,-3.222789e3, 1.575670e4},
5231 {1.189447e4, 5.161815e4, 4.870290e4, 1.266808e4, 2.069603e3}};
5233 for(i = 0; i < sizey; i++) {
5234 for(j = 0; j < sizex; j++) {
5235 egscof[i][j][0] = egs1[i][j];
5236 egscof[i][j][1] = egs2[i][j];
5237 egscof[i][j][2] = egs3[i][j];
5238 egscof[i][j][3] = egs4[i][j];
5243 if (iz < 19 || iz > 111) {
5247 if(iz > 102 && il > 0) {
5254 amin= 1.2e0*
z + 0.01e0*
z*
z;
5255 amax= 5.8e0*z - 0.024e0*z*
z;
5257 if(a < amin || a > amax) {
5269 for(i = 0; i < 7; i++) {
5270 for(j = 0; j < 7; j++) {
5271 bfis0=bfis0+elzcof[j][i]*pz[i]*pa[j];
5283 amin2 = 1.4e0*z + 0.009e0*z*
z;
5284 amax2 = 20.e0 + 3.0e0*
z;
5286 if((a < amin2-5.e0 || a > amax2+10.e0) && il > 0) {
5296 for(i = 0; i < 4; i++) {
5297 for(j = 0; j < 5; j++) {
5298 el80 = el80 + elmcof[i][j]*pz[j]*pa[i];
5299 el20 = el20 + emncof[i][j]*pz[j]*pa[i];
5310 for(i = 0; i < 4; i++) {
5311 for(j = 0; j < 6; j++) {
5312 elmax = elmax + emxcof[i][j]*pz[j]*pa[i];
5323 x = sel20/(*selmax);
5324 y = sel80/(*selmax);
5328 q = 0.2/(std::pow(sel20,2)*std::pow(sel80,2)*(sel20-sel80));
5329 qa = q*(4.0*std::pow(sel80,3) - std::pow(sel20,3));
5330 qb = -q*(4.0*std::pow(sel80,2) - std::pow(sel20,2));
5331 bfis = bfis*(1.0 + qa*std::pow(el,2) + qb*std::pow(el,3));
5335 aj = (-20.0*std::pow(
x,5) + 25.e0*std::pow(
x,4) - 4.0)*std::pow((
y-1.0),2)*
y*
y;
5336 ak = (-20.0*std::pow(y,5) + 25.0*std::pow(y,4) - 1.0) * std::pow((
x-1.0),2)*
x*
x;
5337 q = 0.2/(std::pow((y-x)*((1.0-x)*(1.0-y)*x*y),2));
5338 qa = q*(aj*y - ak*
x);
5339 qb = -q*(aj*(2.0*y + 1.0) - ak*(2.0*x + 1.0));
5341 a1 = 4.0*std::pow(z,5) - 5.0*std::pow(z,4) + 1.0;
5342 a2 = qa*(2.e0*z + 1.e0);
5343 bfis=bfis*(a1 + (z - 1.e0)*(a2 + qb*z)*z*z*(z - 1.e0));
5350 if(el > (*selmax)) {
5356 if(el > (*selmax)) {
5360 for(k = 0; k < 4; k++) {
5361 for(l = 0; l < 6; l++) {
5362 for(
m = 0;
m < 5;
m++) {
5363 egs = egs + egscof[l][
m][k]*pz[l]*pa[k]*pl[2*
m];
5409 ferf=-
gammp(0.5,x*x);
5411 ferf=
gammp(0.5,x*x);;
5421 if(x<0.0 || a<=0.0)std::cout <<
"G4Abla::gammp = bad arguments in gammp" << std::endl;
5423 gser(&gamser,a,x,gln);
5426 gcf(&gammcf,a,x,gln);
5445 for(
G4int i=1;i<=itmax;i++){
5449 if(std::fabs(d)<fpmin)d=fpmin;
5451 if(std::fabs(c)<fpmin)c=fpmin;
5455 if(std::fabs(del-1.)<
eps)
goto dir1;
5457 std::cout <<
"a too large, ITMAX too small in gcf" << std::endl;
5459 fgammcf=std::exp(-x+a*std::log(x)-gln)*h;
5472 if(x<0.)std::cout <<
"G4Abla::gser = x < 0 in gser" << std::endl;
5483 if(std::fabs(del)<std::fabs(sum)*eps)
goto dir1;
5485 std::cout <<
"a too large, ITMAX too small in gser" << std::endl;
5487 fgamser=sum*std::exp(-x+a*std::log(x)-gln);
5495 G4double cof[6]={76.18009172947146,-86.50532032941677,24.01409824083091,
5496 -1.231739572450155,0.1208650973866179e-2,-0.5395239384953e-5};
5502 tmp=(x+0.5)*std::log(tmp)-
tmp;
5503 ser=1.000000000190015;
5504 for(
G4int j=0;j<6;j++){
5509 return fgammln=tmp+std::log(stp*ser/x);
5517 return (E*std::exp(-E));
5523 return (1.0 - (E + 1.0) * std::exp(-E));
5539 const G4int pSize = 101;
5557 for(i = 1; i <= 99; i++) {
5561 if (std::fabs(
f(x) -
G4double(i)/100.0) < 1
e-5) {
5586 x = (p[i] - p[i-1])*(y*100 - i) + p[i];
5603 if(ii <= 0 || jj < 0) {
5614 fpace2=fpace2/1000.;
5616 if(
pace->
dm[ii][jj] == 0.) {
5621 guet(&a, &z, &fpace2);
5622 fpace2=fpace2-ii*931.5;
5623 fpace2=fpace2/1000.;
5642 const G4int qrows = 50;
5643 const G4int qcols = 70;
5645 for(
G4int init_i = 0; init_i < qrows; init_i++) {
5646 for(
G4int init_j = 0; init_j < qcols; init_j++) {
5647 q[init_i][init_j] = 0.0;
5688 G4double aux=1.+(9.*xjj/4./qq/x13);
5690 G4double ee4=avol*xx+asur*(std::pow(xx,.666))+ac*x13+azer;
5691 G4double tota = ee1 + ee2 + ee3 + ee4;
5692 find = 939.55*xneu+938.77*zz - tota;
5704 const G4double fmp = 938.27231,fmn=939.56563;
5710 for(
G4int i=0;i<IMULTBU;i++){
5720 G4double gamma = std::sqrt(1.0 - v2 / (c*c));
5731 for(
G4int i=0;i<IEV_TAB;i++){
5742 G4double gamma = std::sqrt(1.0 - v2 / (c*c));
5810 return -1.0*std::abs(a);
5822 return -1*std::abs(a);
5831 fractpart = std::modf(number, &intpart);
5836 if(fractpart < 0.5) {
5837 return G4int(std::floor(number));
5840 return G4int(std::ceil(number));
5844 if(fractpart < -0.5) {
5845 return G4int(std::floor(number));
5848 return G4int(std::ceil(number));
5852 return G4int(std::floor(number));
5861 mylocaltime = localtime(&mytime);
5864 return(mylocaltime->tm_hour*60*60 + mylocaltime->tm_min*60 + mylocaltime->tm_sec);
5885 if(x-std::floor(x) <= std::ceil(x)-x)
5886 value = double(std::floor(x));
5888 value = double(std::ceil(x));
5897 if(x-std::floor(x) <= std::ceil(x)-x)
5898 value =
G4int(std::floor(x));
5900 value =
G4int(std::ceil(x));
5907 if(x-std::floor(x) <= std::ceil(x)-x)
5908 return G4int(std::floor(x));
5910 return G4int(std::ceil(x));
5915 if(a < b && a < c) {
5918 if(b < a && b < c) {
5921 if(c < a && c < b) {
5941 G4int IZPART,IAPART,NMOTHER;
5944 G4double INT1,INT2,INT3,AKONST,EARG,R0,MPARTNER;
5947 G4double PAR_A1=0.,PAR_B1=0.,FACT=1.;
5956 NMOTHER =
idnint(AMOTHER-ZMOTHER);
5965 HBAR = 6.582122e-22;
5970 APARTNER = AMOTHER - APART;
5971 MPARTNER = APARTNER * 931.49 /
C2;
5974 if(IAPART==1&&IZPART==0){
5976 MPART = 939.56 /
C2;
5978 if(IAPART==1&&IZPART==1){
5980 MPART = 938.27 /
C2;
5983 if(IAPART==2&&IZPART==1){
5985 MPART = 1876.10 /
C2;
5987 if(IAPART==3&&IZPART==1){
5989 MPART = 2809.39 /
C2;
5991 if(IAPART==3&&IZPART==2){
5993 MPART = 2809.37 /
C2;
5995 if(IAPART==4&&IZPART==2){
5997 MPART = 3728.35 /
C2;
6001 MPART = APART * 931.49 /
C2;
6010 MU = MPARTNER * MPART / (MPARTNER + MPART);
6014 RGEOM = R0 * (std::pow(APART,1.0/3.0)+std::pow(AMOTHER-APART,1.0/3.0));
6017 AKONST = HBAR*std::sqrt(1.0 / MU);
6020 BKONST = MPART / ( PI * PI * HBAR * HBAR);
6024 INT1 = 2.0 * std::pow(TEMP,3.) / (2.0 * TEMP +
B);
6026 ARG = std::sqrt(B/TEMP);
6027 EARG = (
erf(ARG) - 1.0);
6028 if(std::abs(EARG)<1.e-9) EARG = 0.0;
6030 INT2 = 0.5 * std::sqrt(PI) * std::pow(TEMP,3.0/2.0);
6033 if(AEXP>700.0) AEXP = 700.0;
6034 INT2 = (2.0*B*B +TEMP*
B)/std::sqrt(B) + std::exp(AEXP) * std::sqrt(PI/(4.0*TEMP))*(4.0*B*B+4.0*B*TEMP - TEMP*TEMP) *EARG;
6035 if(INT2<0.0) INT2 = 0.0;
6038 if(EARG==0.0) INT2 = 0.0;
6041 INT3 = 2.0*TEMP*TEMP*TEMP / (2.0*TEMP*TEMP + 4.0*B*TEMP + B*
B);
6043 if(IZPART<-1.0&&ZMOTHER<151.0){
6047 fwidth = PI * BKONST * G * std::sqrt((RGEOM * RGEOM * INT1 + 2.0 * AKONST * RGEOM * INT2 + AKONST * AKONST * INT3) * RGEOM * RGEOM * INT1);
6050 fwidth = PI * BKONST * G *(RGEOM * RGEOM * INT1 + 2.0 * AKONST * RGEOM * INT2 + AKONST * AKONST * INT3);
6058 PAR_A1=std::exp(2.302585*0.2083*std::exp(-0.01548472*AMOTHER))-0.05;
6059 PAR_B1 = 0.59939389 + 0.00915657 * AMOTHER;
6061 if(AMOTHER>154.0&&AMOTHER<195.0){
6062 PAR_A1=1.0086961-8.629e-5*AMOTHER;
6063 PAR_B1 = 1.5329331 + 0.00302074 * AMOTHER;
6065 if(AMOTHER>194.0&&AMOTHER<208.0){
6066 PAR_A1=9.8356347-0.09294663*AMOTHER+2.441e-4*AMOTHER*AMOTHER;
6067 PAR_B1 = 7.7701987 - 0.02897401 * AMOTHER;
6069 if(AMOTHER>207.0&&AMOTHER<228.0){
6070 PAR_A1=15.107385-0.12414415*AMOTHER+2.7222e-4*AMOTHER*AMOTHER;
6071 PAR_B1=-64.078009+0.56813179*AMOTHER-0.00121078*AMOTHER*AMOTHER;
6074 if(
mod(NMOTHER,2)==0&&NMOTHER>147.){
6075 PAR_A1 = 2.0*(0.9389118 + 6.4559e-5 * AMOTHER);
6077 if(
mod(NMOTHER,2)==1)PAR_A1 = 3.0*(0.9389118 + 6.4559
e-5 * AMOTHER);
6079 PAR_B1 = 2.1507177 + 0.00146119 * AMOTHER;
6085 FACT = std::exp((2.302585*PAR_A1*std::exp(-PAR_B1*(EXC-SB))));
6086 if(FACT<1.0) FACT = 1.0;
6087 if(IZPART<-1.&&ZMOTHER<151.0){
6089 fwidth = fwidth / std::sqrt(FACT);
6091 fwidth = fwidth / FACT;
6096 std::cout <<
"LOOK IN PARTICLE_WIDTH!" << std::endl;
6097 std::cout <<
"ACN,APART :"<< AMOTHER << APART << std::endl;
6098 std::cout <<
"EXC,TEMP,B,SB :" << EXC <<
" " << TEMP <<
" " << B <<
" " << SB << std::endl;
6099 std::cout <<
"INTi, i=1-3 :" << INT1 <<
" " << INT2 <<
" " << INT3 << std::endl;
6100 std::cout <<
" " << std::endl;
6117 MU = (A -
ap) * ap / A;
6121 HO = 197.3287 * omega;
6126 fpen=std::pow(10.0,4.
e-4*std::pow(T/(HO*HO*std::pow(MU,0.25)),-4.3/2.3026));
6144 (*BS) = 1.0 + 0.4*ALPHA2*ALPHA2 - 4.0/105.0*ALPHA2*ALPHA2*ALPHA2 - 66.0/175.0*ALPHA2*ALPHA2*ALPHA2*ALPHA2 - 4.0/35.0*ALPHA2*ALPHA2*ALPHA4 + ALPHA4*ALPHA4;
6146 (*BK) = 1.0 + 0.4*ALPHA2*ALPHA2 + 16.0/105.0*ALPHA2*ALPHA2*ALPHA2 - 82.0/175.0*ALPHA2*ALPHA2*ALPHA2*ALPHA2 + 2.0/35.0*ALPHA2*ALPHA2*ALPHA4 + ALPHA4*ALPHA4;
6190 G4double DEFO_INIT,OMEGA,HOMEGA,OMEGA_GS,HOMEGA_GS,K1,MFCD;
6191 G4double BET1,XACT,SIGMA_SQR,W_EXP,XB,NORM,SIGMA_SQR_INF,W_INFIN,W;
6192 G4double FUNC_TRANS,LOG_SLOPE_INF,LOG_SLOPE_ABS;
6199 fomega_gs(AF,ZF,&K1,&OMEGA_GS,&HOMEGA_GS);
6203 if((bet*bet)>4.0*OMEGA_GS*OMEGA_GS){
6204 BET1=std::sqrt(bet*bet-4.0*OMEGA_GS*OMEGA_GS);
6209 SIGMA_SQR = (FT/K1)*(1.0 -((2.0*bet*bet/(BET1*BET1)* (0.5 * (std::exp(0.50*(BET1-bet)*1.e21*TIME) - std::exp(0.5*(-BET1-bet)*1.e21*TIME)))*(0.5 * (std::exp(0.50*(BET1-bet)*1.e21*TIME) - std::exp(0.5*(-BET1-bet)*1.e21*TIME)))) + (bet/BET1*0.50 * (std::exp((BET1-bet)*1.e21*TIME)-std::exp((-BET1-bet)*1.e21*TIME))) + 1. * std::exp(-bet*1.e21*TIME)));
6212 XACT = DEFO_INIT *std::exp(-0.5*(bet-BET1)*1.e21*(TIME-T_0));
6217 BET1=std::sqrt(4.0*OMEGA_GS*OMEGA_GS-bet*bet);
6218 SIGMA_SQR = FT/K1*(1.-std::exp(-1.0*bet*1.e21*TIME)*(bet*bet/(BET1*BET1)*(1.-std::cos(BET1*1.e21*TIME)) + bet/BET1*std::sin(BET1*1.e21*TIME) + 1.0));
6219 XACT = DEFO_INIT*std::cos(0.5*BET1*1.e21*(TIME-T_0))*std::exp(-bet*1.e21*(TIME-T_0));
6225 XB = 7./3.*Y-938./765.*Y*Y+9.499768*Y*Y*Y-8.050944*Y*Y*Y*
Y;
6230 NORM = 1./std::sqrt(2.*PI*SIGMA_SQR);
6232 W_EXP = -1.*(XB - XACT)*(XB - XACT)/(2.0 * SIGMA_SQR);
6233 if(W_EXP<(-708.0) ) W_EXP = -708.0;
6234 W = NORM * std::exp( W_EXP ) * FT / (K1 * SIGMA_SQR);
6241 SIGMA_SQR_INF = FT/K1;
6242 W_EXP = -XB*XB/(2.0 * SIGMA_SQR_INF);
6243 if(W_EXP<(-708.0))W_EXP = -708.0;
6244 W_INFIN = std::exp(W_EXP)/std::sqrt(2.0*PI*SIGMA_SQR_INF);
6245 FUNC_TRANS = W / W_INFIN;
6250 LOG_SLOPE_INF =
cram(bet,HOMEGA)*bet*MFCD*OMEGA/FT;
6251 LOG_SLOPE_ABS = (XB-XACT)/SIGMA_SQR-XB/SIGMA_SQR_INF+
cram(bet,HOMEGA)*bet*MFCD*OMEGA/FT;
6253 FUNC_TRANS = FUNC_TRANS * LOG_SLOPE_ABS/LOG_SLOPE_INF;
6259 void G4Abla::part_fiss(
G4double BET,
G4double GP,
G4double GF,
G4double Y,
G4double TAUF,
G4double TS1,
G4double TSUM,
G4int *CHOICE,
G4double ZF,
G4double AF,
G4double FT,
G4double *T_LAPSE,
G4double *GF_LOC)
6318 G4double K1,OMEGA,HOMEGA,t_0,STEP_LENGTH,LOC_TIME_BEGIN,LOC_TIME_END=0.,BEGIN_TIME=0.,FISS_PROB,
X,TS2,LAMBDA,REAC_PROB;
6336 if(BET*BET>4.0*OMEGA*OMEGA){
6343 t_0 = BET*1.e21*HBAR*HBAR/(4.*HOMEGA*FT)/16.;
6346 if(((2.*FT-HOMEGA/16.)>0.000001) && BET>0.0){
6351 t_0 = (std::log(2.*FT/(2.*FT-HOMEGA/16.)))/(BET*1.e21);
6363 STEP_LENGTH = 1.5*TAUF/50.;
6368 BEGIN_TIME = TSUM + t_0;
6370 if(BEGIN_TIME<0.0) std::cout <<
"CURRENT TIME < 0" << BEGIN_TIME << std::endl;
6372 if(BEGIN_TIME<1.50*TAUF){
6373 LOC_TIME_BEGIN = BEGIN_TIME;
6375 while((LOC_TIME_BEGIN<1.5*TAUF)&&fchoice==0){
6377 LOC_TIME_END = LOC_TIME_BEGIN + STEP_LENGTH;
6380 fGF_LOC=(
func_trans(LOC_TIME_BEGIN,ZF,AF,BET,Y,FT,t_0)+
func_trans(LOC_TIME_END,ZF,AF,BET,Y,FT,t_0))/2.0;
6382 fGF_LOC = fGF_LOC * GF;
6392 LAMBDA = 1.0/TS1 + 1.0/TS2;
6398 REAC_PROB = std::exp(-1.0*STEP_LENGTH*LAMBDA);
6403 FISS_PROB = fGF_LOC / (fGF_LOC+GP);
6414 LOC_TIME_BEGIN = LOC_TIME_END;
6417 fT_LAPSE = LOC_TIME_END - BEGIN_TIME;
6424 FISS_PROB = GF / (GF+GP);
6434 LAMBDA = 1./TS1 + 1./TS2;
6456 (*T_LAPSE)=fT_LAPSE;
6469 G4double MFCD,OMEGA,HOMEGA1,HOMEGA2=0.,GFTUN;
6470 G4double E1,E2,EXP_FACT,CORR_FUNCT,FACT1,FACT2,FACT3;
6478 if(
mod(IN,2)==0&&
mod(IZ,2)==0){
6481 EE = EE - 12.0/std::sqrt(A);
6485 if(
mod(IN,2)==1&&
mod(IZ,2)==1){
6489 if(
mod(IN,2)==1&&
mod(IZ,2)==0){
6493 if(
mod(IN,2)==0&&
mod(IZ,2)==1){
6497 E1 = EF + HOMEGA1/2.0/PI*std::log(HOMEGA1*(2.0*PI+HOMEGA2)/4.0/PI/PI);
6499 E2 = EF + HOMEGA2/(2.0*PI)*std::log(1.0+2.0*PI/HOMEGA2);
6506 EXP_FACT = (EE-EF)/(HOMEGA2/(2.0*PI));
6507 if(EXP_FACT>700.0) EXP_FACT = 700.0;
6508 CORR_FUNCT = HOMEGA1 * (1.0-1.0/(1.0+std::exp(EXP_FACT)));
6509 if(
mod(IN,2)==0&&
mod(IZ,2)==0){
6510 CORR_FUNCT = HOMEGA1 * (1.0-1.0/(1.0+std::exp(EXP_FACT)));
6513 FACT1 = HOMEGA1/(2.0*PI*TEMP+HOMEGA1);
6514 FACT2 = (2.0*PI/(2.0*PI+HOMEGA2)-HOMEGA1*(2.0*PI+HOMEGA2)/4.0/PI/PI)/(E2-E1);
6515 FACT3 = HOMEGA2/(2.0*PI*TEMP-HOMEGA2);
6518 GFTUN = FACT1*(std::exp(EE/TEMP)*std::exp(2.0*PI*(EE-EF)/HOMEGA1)-std::exp(-2.0*PI*EF/HOMEGA1));
6521 GFTUN = std::exp(EE/TEMP)*(0.50+FACT2*(EE-EF-TEMP))-std::exp(E1/TEMP)*(0.5+FACT2*(E1-EF-TEMP))+FACT1*(std::exp(E1/TEMP)*std::exp(2.0*PI*(E1-EF)/HOMEGA1)-std::exp(-2.0*PI*EF/HOMEGA1));
6523 GFTUN = std::exp(EE/TEMP)*(1.0+FACT3*std::exp(-2.0*PI*(EE-EF)/HOMEGA2))-std::exp(E2/TEMP)*(1.0+FACT3*std::exp(-2.0*PI*(E2-EF)/HOMEGA2))+std::exp(E2/TEMP)*(0.5+FACT2*(E2-EF-TEMP))-std::exp(E1/TEMP)*(0.5+FACT2*(E1-EF-TEMP))+FACT1*(std::exp(E1/TEMP)*std::exp(2.0*PI*(E1-EF)/HOMEGA1)-std::exp(-2.0*PI*EF/HOMEGA1));
6526 GFTUN = GFTUN/std::exp(EE/TEMP)*DENSF*ENH_FACT/DENSG/2.0/PI;
6527 GFTUN = GFTUN * CORR_FUNCT;
6532 void G4Abla::fission_width(
G4double ZPRF,
G4double A,
G4double EE,
G4double BS,
G4double BK,
G4double EF,
G4double Y,
G4double *GF,
G4double *TEMP,
G4double JPR,
G4int IEROT,
G4int FF_ALLOWED,
G4int OPTCOL,
G4int OPTSHP,
G4double DENSG)
6535 G4double FNORM,MASS_ASYM_SADD_B,FP_PER,FP_PAR,SIG_PER_SP,SIG_PAR_SP;
6536 G4double Z2OVERA,ftemp,fgf,DENSF,ECOR,EROT,qr;
6537 G4double DCR,UCR,ENH_FACTA,ENH_FACTB,ENH_FACT,PONFE;
6542 Z2OVERA = ZPRF * ZPRF /
A;
6545 if((ZPRF<=55.0) || (FF_ALLOWED==0)){
6555 densniv(A,ZPRF,EE,EF,&DENSF,0.0,BS,BK,&ftemp,OPTSHP,0,Y,&ECOR,JPR,1,&qr);
6558 fgf= DENSF/DENSG/PI/2.0*ftemp;
6570 FNORM = 1.2*1.2 * 931.49 * 1.e-2 / (9.0 * 6.582122*6.582122);
6573 FP_PER = 2.0/5.0*std::pow(A,5.0/3.0)*FNORM*(1. + 7.0/6.0*Y*(1.0+1396.0/255.*
Y));
6578 if(Z2OVERA<=30.0) FP_PER = 6.50;
6581 FP_PAR = 2.0/5.0*std::pow(A,5.0/3.0)*FNORM*(1.0 - 7.0/3.0*Y*(1.0-389.0/255.0*
Y));
6582 if(FP_PAR<0.0) FP_PAR = 0.0;
6584 EROT = JPR * JPR / (2.0 * std::sqrt(FP_PAR*FP_PAR + FP_PER*FP_PER));
6585 if(IEROT==1) EROT = 0.0;
6588 SIG_PER_SP = std::sqrt(FP_PER * ftemp);
6590 if(SIG_PER_SP<1.0) SIG_PER_SP = 1.0;
6593 SIG_PAR_SP = std::sqrt(FP_PAR * ftemp);
6597 MASS_ASYM_SADD_B = 2.0;
6599 MASS_ASYM_SADD_B = 1.0;
6603 if(Z2OVERA>35.&&Z2OVERA<=(110.*110./298.0)){
6605 ENH_FACTA = std::sqrt(8.0*PI) * SIG_PER_SP*SIG_PER_SP * SIG_PAR_SP;
6607 ENH_FACTB = MASS_ASYM_SADD_B * SIG_PER_SP*SIG_PER_SP;
6609 ENH_FACT = ENH_FACTA * ENH_FACTB / (ENH_FACTA + ENH_FACTB);
6613 ENH_FACT = MASS_ASYM_SADD_B*SIG_PER_SP*SIG_PER_SP;
6616 ENH_FACT = std::sqrt(8.0*PI) * SIG_PER_SP*SIG_PER_SP* SIG_PAR_SP;
6621 PONFE = (ECOR-UCR-EROT)/DCR;
6622 if(PONFE>700.) PONFE = 700.0;
6624 ENH_FACT = 1.0/(1.0+std::exp(PONFE))*ENH_FACT+1.0;
6626 if(ENH_FACT<1.0)ENH_FACT = 1.0;
6627 fgf= DENSF/DENSG/PI/2.0*ftemp*ENH_FACT;
6631 fgf=
tunnelling(A,ZPRF,Y,EE,EF,ftemp,DENSG,DENSF,ENH_FACT);
6643 G4double AFRAGMENT,S4FINAL,ALEVDENS;
6644 G4double THETA_MOTHER,THETA_ORBITAL;
6659 if (EEFINAL<=0.01) EEFINAL = 0.01;
6660 AFRAGMENT = AMOTHER - ADAUGHTER;
6661 ALEVDENS = 0.073*AMOTHER + 0.095*std::pow(AMOTHER,2.0/3.0);
6662 S4FINAL = ALEVDENS * EEFINAL;
6663 if(S4FINAL <= 0.0 || S4FINAL > 100000.){
6664 std::cout<<
"S4FINAL:" << S4FINAL << ALEVDENS << EEFINAL <<
idnint(AMOTHER) <<
idnint(AFRAGMENT) << std::endl;
6666 THETA_MOTHER = 0.0111 * std::pow(AMOTHER,1.66667);
6667 THETA_ORBITAL = 0.0323 / std::pow(AMOTHER,2.) *std::pow(std::pow(AFRAGMENT,0.33333) + std::pow(ADAUGHTER,0.33333),2.) * AFRAGMENT*ADAUGHTER*(AFRAGMENT+ADAUGHTER);
6669 *LORBITAL = -1.* THETA_ORBITAL * (LMOTHER / THETA_MOTHER + std::sqrt(S4FINAL) /(ALEVDENS*LMOTHER));
6671 *SIGMA_LORBITAL = std::sqrt(std::sqrt(S4FINAL) * THETA_ORBITAL / ALEVDENS);
6699 G4int iz=0,
in=0,IZIMF=0,INMI=0,INMA=0,IZCN=0,INCN=0,INIMFMI=0,INIMFMA=0,ILIMMAX=0,INNMAX=0,INMIN=0,IAIMF=0,IZSTOP=3,IZMEM=0,IA=0,INMINMEM=0,INMAXMEM=0,IIA=0;
6700 G4double BS=0,BK=0,BC=0,BSHELL=0,DEFBET=0,DEFBETIMF=0,EROT=0,MAIMF=0,MAZ=0,MARES=0,AIMF_1,OMEGAP=0,fBIMF=0.0,BSIMF=0,A1PAR=0,A2PAR=0,SUM_A,EEDAUG;
6701 G4double DENSCN=0,TEMPCN=0,ECOR=0,IINERT=0,EROTCN=0,WIDTH_IMF=0.0,WIDTH1=0,IMFARG=0,QR=0,QRCN=0,DENSIMF=0,fTIMF=0,fZIMF=0,fAIMF=0.0,NIMF=0,fSBIMF=0;
6702 G4double PI = 3.141592653589793238;
6704 G4double SDWprevious=0,SUMDW_TOT=0,SUM_Z=0,
X=0,SUMDW_N_TOT=0,XX=0;
6712 IZIMFMAX =
idnint(ZCN / 2.0);
6715 std::cout <<
"CHARGE_IMF line 46" << std::endl;
6716 std::cout <<
"Problem: IZIMFMAX < 3 " << std::endl;
6717 std::cout <<
"ZCN,IZIMFMAX," << ZCN <<
"," << IZIMFMAX << std::endl;
6725 bsbkbc(ACN,ZCN,&BS,&BK,&BC);
6727 densniv(ACN,ZCN,EE,0.0,&DENSCN,BSHELL,BS,BK,&TEMPCN,0,0,DEFBET,&ECOR,JPRF,0,&QRCN);
6729 IINERT = 0.4 * 931.49 * 1.16*1.16 * std::pow(ACN,5.0/3.0)*(1.0 + 0.5*std::sqrt(5./(4.*PI))*DEFBET);
6730 EROTCN = JPRF * JPRF * 197.328 * 197.328 /(2. * IINERT);
6732 for(IZIMF=3;IZIMF<=IZIMFMAX;IZIMF++){
6741 INIMFMI =
max(1,INIMFMI-2);
6744 INCN =
idnint(ACN) - IZCN;
6748 INMI =
max(1,INMI-2);
6749 INMIN =
max(INIMFMI,INCN-INMA);
6750 INNMAX =
min(INIMFMA,INCN-INMI);
6752 ILIMMAX =
max(INNMAX,INMIN);
6755 for(
G4int INIMF=INMIN;INIMF<=ILIMMAX;INIMF++){
6756 IAIMF = IZIMF + INIMF;
6757 DW[IZIMF][IAIMF] = 0.0;
6758 AIMF_1 = 1.0*(IAIMF);
6761 mglms(ACN-AIMF_1,ZCN-ZIMF_1,OPTSHPIMF,&MARES);
6762 mglms(AIMF_1,ZIMF_1,OPTSHPIMF,&MAIMF);
6763 mglms(ACN,ZCN,OPTSHPIMF,&MAZ);
6767 SSBIMF[IZIMF][IAIMF] = 1.e37;
6770 SSBIMF[IZIMF][IAIMF] = MAIMF + MARES - MAZ + fBIMF;
6771 BBIMF[IZIMF][IAIMF] = fBIMF;
6777 IINERT = 0.40 * 931.490 * 1.160*1.160 * std::pow(ACN,5.0/3.0)*(std::pow(AIMF_1,5.0/3.0) + std::pow(ACN - AIMF_1,5.0/3.0)) + 931.490 * 1.160*1.160 * AIMF_1 * (ACN-AIMF_1) / ACN *(std::pow(AIMF_1,1.0/3.0) + std::pow(ACN - AIMF_1,1.0/3.0))*(std::pow(AIMF_1,1.0/3.0) + std::pow(ACN - AIMF_1,1.0/3.0));
6779 EROT = JPRF * JPRF * 197.328 * 197.328 /(2.0 * IINERT);
6782 if (EE<(SSBIMF[IZIMF][IAIMF]+EROT) || DENSCN<=0.0){
6791 densniv(ACN,ZCN,EE,SSBIMF[IZIMF][IAIMF],&DENSIMF,0.0,BSIMF,1.0,&fTIMF,0,0,DEFBETIMF,&ECOR,JPRF,2,&QR);
6792 IMFARG = (SSBIMF[IZIMF][IAIMF]+EROTCN-EROT)/fTIMF;
6793 if(IMFARG>200.0) IMFARG = 200.0;
6795 WIDTH1 =
width(ACN,ZCN,AIMF_1,ZIMF_1,fTIMF,fBIMF,SSBIMF[IZIMF][IAIMF],EE-EROT);
6797 WIDTH_IMF = WIDTH1 * std::exp(-IMFARG) * QR / QRCN;
6800 std::cout <<
"GAMMA_IMF=0 -> LOOK IN GAMMA_IMF CALCULATIONS!" << std::endl;
6801 std::cout <<
"ACN,ZCN,AIMF,ZIMF:" <<
idnint(ACN) <<
"," <<
idnint(ZCN) <<
"," <<
idnint(AIMF_1) <<
"," <<
idnint(ZIMF_1) << std::endl;
6802 std::cout <<
"SSBIMF,TIMF :" << SSBIMF[IZIMF][IAIMF] <<
"," << fTIMF << std::endl;
6803 std::cout <<
"DEXP(-IMFARG) = " << std::exp(-IMFARG) << std::endl;
6804 std::cout <<
"WIDTH1 =" << WIDTH1 << std::endl;
6808 SDW[IZIMF] = SDW[IZIMF] + WIDTH_IMF;
6810 DW[IZIMF][IAIMF] = WIDTH_IMF;
6818 SDWprevious = 1.e20;
6821 for(
G4int III_ZIMF=3;III_ZIMF<=IZIMFMAX;III_ZIMF++){
6823 if(SDW[III_ZIMF]==0.0){
6824 IZSTOP = III_ZIMF - 1;
6828 if(SDW[III_ZIMF]>SDWprevious){
6829 IZSTOP = III_ZIMF - 1;
6832 SDWprevious = SDW[III_ZIMF];
6844 A1PAR = std::log10(SDW[IZSTOP]/SDW[IZSTOP-2])/std::log10((1.0*IZSTOP)/(1.0*IZSTOP-2.0));
6845 A2PAR = std::log10(SDW[IZSTOP]) - A1PAR * std::log10(1.0*(IZSTOP));
6846 if(A2PAR>0.)A2PAR=-1.*A2PAR;
6847 if(A1PAR>0.)A1PAR=-1.*A1PAR;
6851 for(
G4int II_ZIMF = IZSTOP;II_ZIMF<=IZIMFMAX;II_ZIMF++){
6852 SDW[II_ZIMF] = std::pow(10.0,A2PAR) * std::pow(1.0*II_ZIMF,A1PAR);
6853 if(SDW[II_ZIMF]<0.0) SDW[II_ZIMF] = 0.0;
6860 for(
G4int I_ZIMF = 3;I_ZIMF<=IZIMFMAX;I_ZIMF++){
6861 SUMDW_TOT = SUMDW_TOT + SDW[I_ZIMF];
6864 std::cout <<
"*********************" << std::endl;
6865 std::cout <<
"IMF function" << std::endl;
6866 std::cout <<
"SUM of decay widths = " << SUMDW_TOT <<
" IZIMFMAX = " << IZIMFMAX << std::endl;
6867 std::cout <<
"IZSTOP = " << IZSTOP << std::endl;
6875 X =
haz(1)*SUMDW_TOT;
6882 for(
G4int IZ = 3;IZ<=IZIMFMAX;IZ++){
6883 SUM_Z = SUM_Z + SDW[IZ];
6896 INMINMEM =
max(1,INMINMEM-2);
6899 INMI =
max(1,INMI-2);
6902 INMINMEM =
max(INMINMEM,INCN-INMA);
6903 INMAXMEM =
min(INMAXMEM,INCN-INMI);
6905 INMAXMEM =
max(INMINMEM,INMAXMEM);
6909 for(
G4int IIINIMF = INMINMEM;IIINIMF<=INMAXMEM;IIINIMF++){
6910 IA = IZMEM + IIINIMF;
6911 if(IZMEM>=3&&IZMEM<=95&&IA>=4&&IA<=250){
6912 SUMDW_N_TOT = SUMDW_N_TOT + DW[IZMEM][IA];
6914 std::cout <<
"CHARGE IMF OUT OF RANGE" << IZMEM <<
", " << IA <<
", " <<
idnint(ACN) <<
", " <<
idnint(ZCN) <<
", " << TEMP << std::endl;
6918 XX =
haz(1)*SUMDW_N_TOT;
6921 for(
G4int IINIMF = INMINMEM;IINIMF<=INMAXMEM; IINIMF++){
6922 IIA = IZMEM + IINIMF;
6924 SUM_A = SUM_A + DW[IZMEM][IIA];
6933 NIMF = fAIMF - fZIMF;
6935 if((ACN-ZCN-NIMF)<=0.0 || (ZCN-fZIMF) <= 0.0){
6936 std::cout <<
"IMF Partner unstable:" << std::endl;
6937 std::cout <<
"System: Acn,Zcn,NCN:" << std::endl;
6938 std::cout <<
idnint(ACN) <<
", " <<
idnint(ZCN) <<
", " <<
idnint(ACN-ZCN) << std::endl;
6939 std::cout <<
"IMF: A,Z,N:" << std::endl;
6940 std::cout <<
idnint(fAIMF) <<
", " <<
idnint(fZIMF) <<
", " <<
idnint(fAIMF-fZIMF) << std::endl;
6941 std::cout <<
"Partner: A,Z,N:" << std::endl;
6942 std::cout <<
idnint(ACN-fAIMF) <<
", " <<
idnint(ZCN-fZIMF) <<
", " <<
idnint(ACN-ZCN-NIMF) << std::endl;
6943 std::cout <<
"----nmin,nmax" << INMINMEM <<
", " << INMAXMEM << std::endl;
6944 std::cout <<
"----- warning: Zimf=" << fZIMF <<
" Aimf=" << fAIMF << std::endl;
6945 std::cout <<
"----- look in subroutine IMF" << std::endl;
6946 std::cout <<
"ACN,ZCN,ZIMF,AIMF,temp,EE,JPRF::" << ACN <<
", " << ZCN <<
", " << fZIMF <<
", " << fAIMF <<
", " << TEMP <<
", " << EE <<
", " << JPRF << std::endl;
6947 std::cout <<
"-IZSTOP,IZIMFMAX:" << IZSTOP <<
", " << IZIMFMAX << std::endl;
6948 std::cout <<
"----X,SUM_Z,SUMDW_TOT:" <<
X <<
", " << SUM_Z <<
", " << SUMDW_TOT << std::endl;
6952 if(fZIMF>=ZCN || fAIMF>=ACN || fZIMF<=2 || fAIMF<=3){
6953 std::cout <<
"----nmin,nmax" << INMINMEM <<
", " << INMAXMEM << std::endl;
6954 std::cout <<
"----- warning: Zimf=" << fZIMF <<
" Aimf=" << fAIMF << std::endl;
6955 std::cout <<
"----- look in subroutine IMF" << std::endl;
6956 std::cout <<
"ACN,ZCN,ZIMF,AIMF,temp,EE,JPRF:" << ACN <<
", " << ZCN <<
", " << fZIMF <<
", " << fAIMF <<
", " << TEMP <<
", " << EE <<
", " << JPRF << std::endl;
6957 std::cout <<
"-IZSTOP,IZIMFMAX:" << IZSTOP <<
", " << IZIMFMAX << std::endl;
6958 std::cout <<
"----X,SUM_Z,SUMDW_TOT:" <<
X <<
", " << SUM_Z <<
", " << SUMDW_TOT << std::endl;
6959 for(
G4int III_ZIMF=3;III_ZIMF<=IZIMFMAX;III_ZIMF++)
6960 std::cout <<
"-**Z,SDW:" << III_ZIMF <<
", " << SDW[III_ZIMF] << std::endl;
6970 if((ZCN-fZIMF)<=0.0)std::cout <<
"CHARGE_IMF ZIMF > ZCN" << std::endl;
6971 if((ACN-fAIMF)<=0.0)std::cout <<
"CHARGE_IMF AIMF > ACN" << std::endl;
6976 EEDAUG = (EE - fSBIMF) * (ACN - fAIMF) / ACN;
6977 bsbkbc(ACN - fAIMF,ZCN-fZIMF,&BS,&BK,&BC);
6978 densniv(ACN-fAIMF,ZCN-fZIMF,EEDAUG,0.0,&DENSIMF,BSHELL,BS,BK,&fTIMF,0,0,DEFBET,&ECOR,0.0,0,&QR);
6981 std::cout <<
"----- warning: EE=" << EE <<
"," <<
" S+Bimf=" << fSBIMF << std::endl;
6982 std::cout <<
"----- look in subroutine IMF" << std::endl;
6983 std::cout <<
"IMF will be resampled" << std::endl;
6997 G4int VISOSTAB[191][2]={
7108 *nmin = VISOSTAB[z-1][0];
7109 *nmax = VISOSTAB[z-1][1];
7125 G4double epsiln = 0.0, probp = 0.0, probd = 0.0, probt = 0.0, probn = 0.0, probhe = 0.0, proba = 0.0, probg = 0.0, probimf=0.0, ptotl = 0.0, tcn = 0.0;
7126 G4double sn = 0.0, sbp = 0.0, sbd = 0.0, sbt = 0.0, sbhe = 0.0, sba = 0.0,
x = 0.0, amoins = 0.0, zmoins = 0.0,
sp= 0.0,sd= 0.0,st= 0.0,she= 0.0,sa= 0.0;
7127 G4double ecn = 0.0, ecp = 0.0, ecd = 0.0, ect = 0.0,eche = 0.0,eca = 0.0, ecg = 0.0,
bp = 0.0, bd = 0.0, bt = 0.0, bhe = 0.0, ba = 0.0;
7129 G4double xcv=0.,ycv=0.,zcv=0.,VXOUT=0.,VYOUT=0.,VZOUT=0.;
7131 G4double jprfn=0.0, jprfp=0.0, jprfd=0.0, jprft=0.0, jprfhe=0.0, jprfa=0.0;
7132 G4double ctet1 = 0.0, stet1 = 0.0, phi1 = 0.0;
7135 G4int itest = 0, sortie=0;
7141 G4double time,tauf,tau0,
a0,a1,emin,ts1,tsum=0.;
7142 G4int inttype=0,inum=0,gammadecay = 0;
7151 const G4double mu2 = 931.494*931.494;
7171 a0 = 0.66482503 - 3.4678935 * std::exp(-0.0104002*ee);
7172 a1 = 5.6846e-04 + 0.00574515 * std::exp(-0.01114307*ee);
7173 tauf = (a0 + a1 * zf*zf/std::pow(af,0.3333333)) * tau0;
7176 direct(zf,af,ee,0.,&probp,&probd,&probt,&probn,&probhe,&proba,&probg,&probimf,&probf,&ptotl,
7177 &sn,&sbp,&sbd,&sbt,&sbhe,&sba,
7178 &ecn,&ecp,&ecd,&ect,&eche,&eca,&ecg,
7179 &
bp,&bd,&bt,&bhe,&ba,&
sp,&sd,&st,&she,&sa,&ef,&ts1,inttype,inum,itest,&sortie,&tcn,
7180 &jprfn, &jprfp, &jprfd, &jprft, &jprfhe, &jprfa, &tsum);
7184 if(ptotl<=0.)
goto post100;
7188 if(emin>1e30)std::cout <<
"ERROR AT THE EXIT OF EVAPORA,E>1.D30,AF" << std::endl;
7195 pc = std::sqrt(std::pow((1.0 + ecn/9.3956e2),2.) - 1.0) * 9.3956e2;
7198 else if(probp!=0.0){
7202 pc = std::sqrt(std::pow((1.0 + ecp/9.3827e2),2.) - 1.0) * 9.3827e2;
7205 else if(probd!=0.0){
7209 pc = std::sqrt(std::pow((1.0 + ecd/1.875358e3),2) - 1.0) * 1.875358e3;
7212 else if(probt!=0.0){
7216 pc = std::sqrt(std::pow((1.0 + ect/2.80828e3),2) - 1.0) * 2.80828e3;
7219 else if(probhe!=0.0){
7222 epsiln = she + eche;
7223 pc = std::sqrt(std::pow((1.0 + eche/2.80826e3),2) - 1.0) * 2.80826e3;
7226 else{
if(proba!=0.0){
7230 pc = std::sqrt(std::pow((1.0 + eca/3.72834e3),2) - 1.0) * 3.72834e3;
7252 pc = std::sqrt(std::pow((1.0 + eca/3.72834e3),2) - 1.0) * 3.72834e3;
7255 else if (
x < proba+probhe) {
7259 epsiln = she + eche;
7260 pc = std::sqrt(std::pow((1.0 + eche/2.80826e3),2) - 1.0) * 2.80826e3;
7263 else if (
x < proba+probhe+probt) {
7268 pc = std::sqrt(std::pow((1.0 + ect/2.80828e3),2) - 1.0) * 2.80828e3;
7271 else if (
x < proba+probhe+probt+probd) {
7276 pc = std::sqrt(std::pow((1.0 + ecd/1.875358e3),2) - 1.0) * 1.875358e3;
7279 else if (
x < proba+probhe+probt+probd+probp) {
7284 pc = std::sqrt(std::pow((1.0 + ecp/9.3827e2),2) - 1.0) * 9.3827e2;
7287 else if (
x < proba+probhe+probt+probd+probp+probn) {
7292 pc = std::sqrt(std::pow((1.0 + ecn/9.3956e2),2.) - 1.0) * 9.3956e2;
7295 else if (
x < proba+probhe+probt+probd+probp+probn+probg) {
7302 if(probp==0.0 && probn==0.0 && probd==0.0 && probt==0.0 && proba==0.0 && probhe==0.0 && probimf==0.0 && probf==0.0){
7313 if(gammadecay==1 && ee<=0.01+epsiln){
7322 if(ee<=0.01) ee = 0.010;
7324 if(af<2.5)
goto post100;
7332 ctet1 = 2.0*rnd - 1.0;
7333 stet1 = std::sqrt(1.0 - std::pow(ctet1,2));
7335 phi1 = rnd*2.0*3.141592654;
7336 xcv = stet1*std::cos(phi1);
7337 ycv = stet1*std::sin(phi1);
7342 G4double ETOT_LP = std::sqrt(pc*pc + amoins*amoins * mu2);
7355 &VXOUT,&VYOUT,&VZOUT);
7365 G4double gamma = 1.0/std::sqrt(1.0 - v2 / (c*c));
7366 G4double etot_lp = amoins*mu * gamma;
7376 pteva = std::sqrt(pxeva*pxeva + pyeva*pyeva);
7378 etot = std::sqrt ( pleva*pleva + pteva*pteva + af*af * mu2 );
7379 vx_eva = c * pxeva / etot;
7380 vy_eva = c * pyeva / etot;
7381 vz_eva = c * pleva / etot;
7383 IEV_TAB_SSC = IEV_TAB_SSC +1;
7385 if(time<tauf)
goto post10;
7391 *E_scission_post = ee;
7398 void G4Abla::unbound(
G4double SN,
G4double SP,
G4double SD,
G4double ST,
G4double SHE,
G4double SA,
G4double BP,
G4double BD,
G4double BT,
G4double BHE,
G4double BA,
G4double *PROBF,
G4double *PROBN,
G4double *PROBP,
G4double *PROBD,
G4double *PROBT,
G4double *PROBHE,
G4double *PROBA,
G4double *PROBIMF,
G4double *PROBG,
G4double *ECN,
G4double *ECP,
G4double *ECD,
G4double *ECT,
G4double *ECHE,
G4double *ECA)
7407 e =
dmin1(SBHE,SN,e);
7408 e =
dmin1(SBHE,SBA,e);
7429 *ECP = (-1.0)*SP + BP;
7446 *ECD = (-1.0)*SD + BD;
7463 *ECT = (-1.0)*ST + BT;
7480 *ECHE= (-1.0)*SHE + BHE;
7498 *ECA = (-1.0)*SA + BA;
7578 Delta_U1_shell_max = -2.45;
7584 Delta_U2_shell = -2.45;
7594 G4double Fwidth_asymm1,Fwidth_asymm2,Fwidth_symm;
7596 Fwidth_asymm1 = 0.65;
7597 Fwidth_asymm2 = 0.65;
7639 Friction_factor = 1.0;
7642 G4double cN_asymm1_shell, cN_asymm2_shell;
7643 G4double gamma,gamma_heavy1,gamma_heavy2;
7659 G4double Sasymm1=0.,Sasymm2=0.,Ssymm=0.,Ysum=0.,Yasymm=0.;
7661 G4double wNasymm1_saddle, wNasymm2_saddle, wNsymm_saddle;
7662 G4double wNasymm2_scission, wNsymm_scission;
7663 G4double wNasymm1, wNasymm2, wNsymm;
7686 G4double A_levdens_heavy1,A_levdens_heavy2;
7690 G4double epsilon_1_saddle,epsilon0_1_saddle;
7691 G4double epsilon_2_saddle,epsilon0_2_saddle,epsilon_symm_saddle;
7696 G4double E_eff1_saddle,E_eff2_saddle;
7697 G4double Epot0_mode1_saddle,Epot0_mode2_saddle,Epot0_symm_saddle;
7698 G4double Epot_mode1_saddle,Epot_mode2_saddle,Epot_symm_saddle;
7699 G4double E_defo,E_defo1,E_defo2,E_scission_pre,E_scission_post;
7705 G4double MassCurv_scis, MassCurv_sadd;
7707 G4double Nheavy1_shell,Nheavy2_shell;
7712 G4double Z_scission,N_scission,A_scission;
7716 G4double DSN132,Delta_U1_shell,E_eff0_saddle;
7721 cN_asymm1_shell = 0.700 * N/
Z;
7722 cN_asymm2_shell = 0.040 * N/
Z;
7726 DSN132 = Nheavy1_in - N/Z * Zheavy1_in;
7727 Aheavy1 = Nheavy1_in + Zheavy1_in + 0.340 * DSN132;
7736 Delta_U1_shell = Delta_U1_shell_max + U1NZ_SLOPE * std::abs(DSN132);
7737 Delta_U1_shell =
min(0.,Delta_U1_shell);
7741 Nheavy1 = N/A * Aheavy1;
7742 Aheavy2 = Nheavy2 * A/
N;
7746 A_levdens = A / xLevdens;
7747 gamma = A_levdens / (0.40 * std::pow(A,1.3333)) * FGAMMA;
7748 A_levdens_heavy1 = Aheavy1 / xLevdens;
7749 gamma_heavy1 = A_levdens_heavy1 / (0.40 * std::pow(Aheavy1,1.3333)) * FGAMMA * FGAMMA1;
7750 A_levdens_heavy2 = Aheavy2 / xLevdens;
7751 gamma_heavy2 = A_levdens_heavy2 / (0.40 * std::pow(Aheavy2,1.3333)) * FGAMMA;
7755 E_saddle_scission = (-24. + 0.02227 * Z*Z/std::pow(A,0.33333))*Friction_factor;
7756 E_saddle_scission =
max( 0.0, E_saddle_scission );
7762 Z2_over_A_eff = Z*Z/
A;
7764 if( Z2_over_A_eff< 34.0 )
7765 MassCurv_scis = std::pow(10., -1.093364 + 0.082933 * Z2_over_A_eff - 0.0002602 * Z2_over_A_eff*Z2_over_A_eff);
7767 MassCurv_scis = std::pow(10., 3.053536 - 0.056477 * Z2_over_A_eff+ 0.0002454 * Z2_over_A_eff*Z2_over_A_eff );
7772 MassCurv_sadd = X_s2s * MassCurv_scis;
7774 cN_symm = 8.0 / std::pow(N,2.) * MassCurv_scis;
7775 cN_symm_sadd = 8.0 / std::pow(N,2.) * MassCurv_sadd;
7777 Nheavy1_shell = Nheavy1;
7780 Nheavy1_eff = (cN_symm_sadd*Nsymm + cN_asymm1_shell *
7781 Uwash(E/A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1) *
7785 Uwash(E/A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1));
7787 Nheavy1_eff = (cN_symm_sadd*Nsymm +
7788 cN_asymm1_shell*Nheavy1_shell)
7793 Nheavy2_NZ = Nheavy2;
7794 Nheavy2_shell = Nheavy2_NZ;
7796 Nheavy2_eff = (cN_symm_sadd*Nsymm +
7798 Uwash(E/A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2) *
7802 Uwash(E/A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2));
7804 Nheavy2_eff = (cN_symm_sadd*Nsymm +
7805 cN_asymm2_shell*Nheavy2_shell)
7809 Delta_U1 = Delta_U1_shell + (Nheavy1_shell - Nheavy1_eff)*(Nheavy1_shell - Nheavy1_eff) * cN_asymm1_shell;
7810 Delta_U1 =
min(Delta_U1,0.0);
7811 Delta_U2 = Delta_U2_shell + (Nheavy2_shell - Nheavy2_eff)*(Nheavy2_shell - Nheavy2_eff) * cN_asymm2_shell;
7812 Delta_U2 =
min(Delta_U2,0.0);
7816 Epot0_mode1_saddle = (Nheavy1_eff-Nsymm)*(Nheavy1_eff-Nsymm) * cN_symm_sadd;
7817 Epot0_mode2_saddle = (Nheavy2_eff-Nsymm)*(Nheavy2_eff-Nsymm) * cN_symm_sadd;
7818 Epot0_symm_saddle = 0.0;
7822 Epot_mode1_saddle = Epot0_mode1_saddle + Delta_U1;
7823 Epot_mode2_saddle = Epot0_mode2_saddle + Delta_U2;
7824 Epot_symm_saddle = Epot0_symm_saddle;
7827 dUeff =
min( Epot_mode1_saddle, Epot_mode2_saddle);
7828 dUeff =
min( dUeff, Epot_symm_saddle);
7829 dUeff = dUeff - Epot_symm_saddle;
7838 epsilon0_1_saddle = Eld - Epot0_mode1_saddle;
7839 epsilon0_2_saddle = Eld - Epot0_mode2_saddle;
7842 epsilon_1_saddle = Eld - Epot_mode1_saddle;
7843 epsilon_2_saddle = Eld - Epot_mode2_saddle;
7845 epsilon_symm_saddle = Eld - Epot_symm_saddle;
7848 eexc1_saddle = epsilon_1_saddle;
7849 eexc2_saddle = epsilon_2_saddle;
7852 EEXC_MAX =
max( eexc1_saddle, eexc2_saddle);
7853 EEXC_MAX =
max( EEXC_MAX, Eld);
7856 epsilon_1_scission = Eld + E_saddle_scission - Epot_mode1_saddle;
7857 epsilon_2_scission = Eld + E_saddle_scission - Epot_mode2_saddle;
7860 epsilon_symm_scission = Eld + E_saddle_scission - Epot_symm_saddle;
7863 E_eff1_saddle = epsilon0_1_saddle - Delta_U1 *
7864 Uwash(epsilon_1_saddle/A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1);
7866 if( E_eff1_saddle < A_levdens * hbom1*hbom1)
7867 E_eff1_saddle = A_levdens * hbom1*hbom1;
7870 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_eff1_saddle) /
7872 Uwash(epsilon_1_saddle/A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)+
7875 E_eff2_saddle = epsilon0_2_saddle -
7877 Uwash(epsilon_2_saddle/A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2);
7879 if(E_eff2_saddle < A_levdens * hbom2*hbom2)
7880 E_eff2_saddle = A_levdens * hbom2*hbom2;
7883 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_eff2_saddle) /
7885 Uwash(epsilon_2_saddle/A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)+
7888 E_eff0_saddle = epsilon_symm_saddle;
7889 if(E_eff0_saddle < A_levdens * hbom3*hbom3)
7890 E_eff0_saddle = A_levdens * hbom3*hbom3;
7893 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_eff0_saddle) /
7896 if(epsilon_symm_scission > 0.0 ){
7897 E_HELP =
max(E_saddle_scission,epsilon_symm_scission);
7899 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*(E_HELP)) /
7903 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_saddle_scission) /
7910 if( E_saddle_scission == 0.0 ){
7911 wNasymm1_scission = wNasymm1_saddle;
7912 wNasymm2_scission = wNasymm2_saddle;
7914 if( Nheavy1_eff > 75.0 ){
7915 wNasymm1_scission = std::sqrt(21.0)*N/
A;
7916 wNasymm2_scission =
max( 12.8 - 1.0 *(92.0 - Nheavy2_eff),1.0)*N/
A;
7919 wNasymm1_scission = wNasymm1_saddle;
7920 wNasymm2_scission = wNasymm2_saddle;
7924 wNasymm1_scission =
max( wNasymm1_scission, wNasymm1_saddle );
7925 wNasymm2_scission =
max( wNasymm2_scission, wNasymm2_saddle );
7927 wNasymm1 = wNasymm1_scission * Fwidth_asymm1;
7928 wNasymm2 = wNasymm2_scission * Fwidth_asymm2;
7929 wNsymm = wNsymm_scission * Fwidth_symm;
7932 Aheavy1_eff = Nheavy1_eff * A/
N;
7933 Aheavy2_eff = Nheavy2_eff * A/
N;
7935 A_levdens_heavy1 = Aheavy1_eff / xLevdens;
7936 A_levdens_heavy2 = Aheavy2_eff / xLevdens;
7937 gamma_heavy1 = A_levdens_heavy1 / (0.40 * std::pow(Aheavy1_eff,1.3333)) * FGAMMA * FGAMMA1;
7938 gamma_heavy2 = A_levdens_heavy2 / (0.40 * std::pow(Aheavy2_eff,1.3333)) * FGAMMA;
7940 if( epsilon_symm_saddle < A_levdens * hbom3*hbom3)
7941 Ssymm = 2.0 * std::sqrt(A_levdens*A_levdens * hbom3*hbom3) +
7942 (epsilon_symm_saddle - A_levdens * hbom3*hbom3)/hbom3;
7944 Ssymm = 2.0 * std::sqrt(A_levdens*epsilon_symm_saddle);
7948 if( epsilon0_1_saddle < A_levdens * hbom1*hbom1 )
7949 Ssymm_mode1 = 2.0 * std::sqrt(A_levdens*A_levdens * hbom1*hbom1) +
7950 (epsilon0_1_saddle - A_levdens * hbom1*hbom1)/hbom1;
7952 Ssymm_mode1 = 2.0 * std::sqrt( A_levdens*epsilon0_1_saddle );
7954 if( epsilon0_2_saddle < A_levdens * hbom2*hbom2 )
7955 Ssymm_mode2 = 2.0 * std::sqrt(A_levdens*A_levdens * hbom2*hbom2) +
7956 (epsilon0_2_saddle - A_levdens * hbom2*hbom2)/hbom2;
7958 Ssymm_mode2 = 2.0 * std::sqrt(A_levdens*epsilon0_2_saddle);
7961 if( epsilon0_1_saddle -
7963 Uwash(epsilon_1_saddle/A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)
7964 < A_levdens * hbom1*hbom1 )
7965 Sasymm1 = 2.0 * std::sqrt( A_levdens*A_levdens * hbom1*hbom1 ) +
7966 (epsilon0_1_saddle - Delta_U1 *
7967 Uwash(epsilon_1_saddle/A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)
7968 - A_levdens * hbom1*hbom1)/hbom1;
7970 Sasymm1 = 2.0 *std::sqrt( A_levdens*(epsilon0_1_saddle - Delta_U1 *
7971 Uwash(epsilon_1_saddle/A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)));
7973 if( epsilon0_2_saddle -
7975 Uwash(epsilon_2_saddle/A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)
7976 < A_levdens * hbom2*hbom2 )
7977 Sasymm2 = 2.0 * std::sqrt( A_levdens*A_levdens * hbom2*hbom2 ) +
7978 (epsilon0_1_saddle-Delta_U1 *
7979 Uwash(epsilon_2_saddle/A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)
7980 - A_levdens * hbom2*hbom2)/hbom2;
7983 std::sqrt( A_levdens*(epsilon0_2_saddle - Delta_U2 *
7984 Uwash(epsilon_2_saddle/A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)));
7986 Yasymm1 = ( std::exp(Sasymm1 - Ssymm) - std::exp(Ssymm_mode1 - Ssymm) ) *
7987 wNasymm1_saddle / wNsymm_saddle * 2.0;
7989 Yasymm2 = ( std::exp(Sasymm2 - Ssymm) - std::exp(Ssymm_mode2 - Ssymm) ) *
7990 wNasymm2_saddle / wNsymm_saddle * 2.0;
7992 Ysum = Ysymm + Yasymm1 + Yasymm2;
7995 Ysymm = Ysymm / Ysum;
7996 Yasymm1 = Yasymm1 / Ysum;
7997 Yasymm2 = Yasymm2 / Ysum;
7998 Yasymm = Yasymm1 + Yasymm2;
8004 if( (epsilon_symm_saddle < epsilon_1_saddle) &&
8005 (epsilon_symm_saddle < epsilon_2_saddle) )
8008 if( epsilon_1_saddle < epsilon_2_saddle )
8016 r_e_o = std::pow(10.0,-0.0170 * (E_saddle_scission + Eld)*(E_saddle_scission + Eld));
8034 if( rmode < Yasymm1 )
8037 if( (rmode > Yasymm1) && (rmode < Yasymm) )
8046 N1mean = Nheavy1_eff;
8050 N1mean = Nheavy2_eff;
8066 while( N1r < 5.0 || N2r < 5.0 ){
8085 E_scission_pre =
max( epsilon_1_scission, 1.0 );
8087 if( N1mean > N*0.50 ){
8097 E_scission_pre =
max( epsilon_2_scission, 1.0 );
8098 if( N1mean > N*0.50 ){
8099 beta1 = (N1r - 92.0) * 0.030 + 0.60;
8104 beta1 =
max(beta1,beta1gs);
8105 beta2 = 1.0 - beta1;
8106 beta2 =
max(beta2,beta2gs);
8112 beta2 = (N2r -92.0) * 0.030 + 0.60;
8113 beta2 =
max(beta2,beta2gs);
8114 beta1 = 1.0 - beta2;
8115 beta1 =
max(beta1,beta1gs);
8130 beta =
max(0.177963+0.0153241*Zsymm-1.62037
e-4*Zsymm*Zsymm,betags);
8131 beta1 =
max(0.177963+0.0153241*Z1UCD-1.62037
e-4*Z1UCD*Z1UCD,beta1gs);
8132 beta2 =
max(0.177963+0.0153241*Z2UCD-1.62037
e-4*Z2UCD*Z2UCD,beta2gs);
8134 E_asym =
frldm( Z1UCD, N1r, beta1 ) +
8135 frldm( Z2UCD, N2r, beta2 ) +
8136 ecoul( Z1UCD, N1r, beta1, Z2UCD, N2r, beta2, 2.0 ) -
8137 2.0 *
frldm( Zsymm, Nsymm, beta ) -
8138 ecoul( Zsymm, Nsymm, beta, Zsymm, Nsymm, beta, 2.0 );
8139 E_scission_pre =
max( epsilon_symm_scission - E_asym, 1. );
8148 if(E_scission_pre>5.){
8150 &A_scission,&Z_scission,vx_eva_sc,vy_eva_sc,vz_eva_sc);
8151 N_scission = A_scission - Z_scission;
8155 E_scission_post = E_scission_pre;
8156 N_scission = A_scission - Z_scission;
8162 N1r = N1r * N_scission /
N;
8163 N2r = N2r * N_scission /
N;
8164 Z1UCD = Z1UCD * Z_scission /
Z;
8165 Z2UCD = Z2UCD * Z_scission /
Z;
8203 CZ = (
frldm( Z1UCD-1.0, N1r+1.0, beta1 ) +
8204 frldm( Z2UCD+1.0, N2r-1.0, beta2 ) +
8205 frldm( Z1UCD+1.0, N1r-1.0, beta1 ) +
8206 frldm( Z2UCD-1.0, N2r+1.0, beta2 ) +
8207 ecoul( Z1UCD-1.0, N1r+1.0, beta1,
8208 Z2UCD+1.0, N2r-1.0, beta2, 2.0) +
8209 ecoul( Z1UCD+1.0, N1r-1.0, beta1,
8210 Z2UCD-1.0, N2r+1.0, beta2, 2.0) -
8211 2.0*
ecoul( Z1UCD, N1r, beta1, Z2UCD, N2r, beta2, 2.0) -
8212 2.0*
frldm( Z1UCD, N1r, beta1 ) -
8213 2.0*
frldm( Z2UCD, N2r, beta2) ) * 0.50;
8215 if(1.0/A_levdens*E_scission_post < 0.0)
8216 std::cout <<
"DSQRT 1 < 0" << A_levdens <<
" " << E_scission_post << std::endl;
8218 if(0.50 * std::sqrt(1.0/A_levdens*E_scission_post) / CZ < 0.0){
8219 std::cout <<
"DSQRT 2 < 0 " << CZ << std::endl;
8220 std::cout <<
"This event was not considered" << std::endl;
8224 ZA1width = std::sqrt(0.5*std::sqrt(1.0/A_levdens*E_scission_post)/CZ);
8235 ZA1width =
max(ZA1width,sigZmin);
8237 if(imode == 1 && cpol1 != 0.0){
8241 Z1rr = Z1UCD - cpol1 * A_scission/N_scission;
8247 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED" << std::endl;
8250 if ((
utilabs(Z1rr - Z1r) > 3.0*ZA1width) || Z1r<1.0)
goto fiss2801;
8253 if( imode == 2 && cpol2 != 0.0 ){
8257 Z1rr = Z1UCD - cpol2 * A_scission/N_scission;
8262 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED" << std::endl;
8265 if( (
utilabs(Z1rr - Z1r) > 3.0*ZA1width) || Z1r < 1.0 )
goto fiss2802;
8273 re1 =
frldm( Z1UCD-1.0, N1r+1.0, beta1 ) +
8274 frldm( Z2UCD+1.0, N2r-1.0, beta2 ) +
8275 ecoul( Z1UCD-1.0, N1r+1.0, beta1,
8276 Z2UCD+1.0, N2r-1.0, beta2, d );
8277 re2 =
frldm( Z1UCD, N1r, beta1) +
8278 frldm( Z2UCD, N2r, beta2 ) +
8279 ecoul( Z1UCD, N1r, beta1,
8280 Z2UCD, N2r, beta2, d );
8281 re3 =
frldm( Z1UCD+1.0, N1r-1.0, beta1 ) +
8282 frldm( Z2UCD-1.0, N2r+1.0, beta2 ) +
8283 ecoul( Z1UCD+1.0, N1r-1.0, beta1,
8284 Z2UCD-1.0, N2r+1.0, beta2, d );
8285 eps2 = ( re1 - 2.0*re2 + re3 ) / 2.0;
8286 eps1 = ( re3 - re1 ) / 2.0;
8287 DN1_POL = -eps1 / ( 2.0 * eps2 );
8289 Z1rr = Z1UCD + DN1_POL;
8294 DN1_POL = DN1_POL - 0.6 *
Uwash(E_scission_post,Ecrit,FREDSHELL,gamma);
8295 Z1rr = Z1UCD + DN1_POL;
8296 if ( Z1rr < 50. ) Z1rr = 50.0;
8298 DN1_POL = DN1_POL + 0.60 *
Uwash(E_scission_post,Ecrit,FREDSHELL,gamma);
8299 Z1rr = Z1UCD + DN1_POL;
8300 if ( Z1rr > 50.0 ) Z1rr = 50.0;
8310 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED" << std::endl;
8314 if( (
utilabs(Z1rr - Z1r) > 3.0*ZA1width) || (Z1r < 1.0) )
goto fiss2803;
8326 z2 =
dint( Z_scission ) - z1;
8328 N2 =
dint( N_scission ) - N1;
8332 if( (z1 < 0) || (z2 < 0) || (a1 < 0) || (a2 < 0) ){
8333 std::cout <<
" -------------------------------" << std::endl;
8334 std::cout <<
" Z, A, N : " << Z <<
" " << A <<
" " << N << std::endl;
8335 std::cout << z1 <<
" " << z2 <<
" " << a1 <<
" " << a2 << std::endl;
8336 std::cout << E_scission_post <<
" " << A_levdens <<
" " << CZ << std::endl;
8338 std::cout <<
" -------------------------------" << std::endl;
8347 if( N1mean > N*0.50 ){
8351 if(beta2< beta2gs) beta2 = beta2gs;
8352 E1exc = E_scission_pre * a1 / A + E_defo;
8353 E_defo =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8354 E2exc = E_scission_pre * a2 / A + E_defo;
8358 if(beta1< beta1gs) beta1 = beta1gs;
8359 E_defo =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8360 E1exc = E_scission_pre * a1 / A + E_defo;
8362 E2exc = E_scission_pre * a2 / A + E_defo;
8369 if( N1mean > N*0.5 ){
8372 if(beta1< beta1gs) beta1 = beta1gs;
8373 E_defo =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8374 E1exc = E_scission_pre * a1 / A + E_defo;
8376 if(beta2< beta2gs) beta2 = beta2gs;
8377 E_defo =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8378 E2exc = E_scission_pre * a2 / A + E_defo;
8382 if(beta2< beta2gs) beta2 = beta2gs;
8383 E_defo =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8384 E2exc = E_scission_pre * a2 / A + E_defo;
8386 if(beta1< beta1gs) beta1 = beta1gs;
8387 E_defo =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8388 E1exc = E_scission_pre * a1 / A + E_defo;
8395 if(beta1< beta1gs) beta1 = beta1gs;
8397 if(beta2< beta2gs) beta2 = beta2gs;
8398 E_defo1 =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8399 E_defo2 =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8400 E1exc = E_scission_pre * a1 / A + E_defo1;
8401 E2exc = E_scission_pre * a2 / A + E_defo2;
8406 TKER = ( z1 * z2 * 1.440 ) /
8407 ( R0 * std::pow(a1,0.333330) * (1.0 + 2.0/3.0 * beta1 ) +
8408 R0 * std::pow(a2,0.333330) * (1.0 + 2.0/3.0 * beta2 ) + 2.0 );
8410 EkinR1 = TKER * a2 /
A;
8411 EkinR2 = TKER * a1 /
A;
8412 v1 = std::sqrt(EkinR1/a1) * 1.3887;
8413 v2 = std::sqrt(EkinR2/a2) * 1.3887;
8421 e1 =
gausshaz(0,E1exc,E1exc_sigma);
8422 if(e1<0.)
goto fis987;
8425 e2 =
gausshaz(0,E2exc,E2exc_sigma);
8426 if(e2<0.)
goto fis988;
8465 G4double r_in = 0.0, r_rest = 0.0, r_help = 0.0;
8471 r_in = r_origin + 0.5;
8473 if (r_even_odd < 0.001) {
8474 i_out = (
G4int)(r_floor);
8477 r_rest = r_in - r_floor;
8478 r_middle = r_floor + 0.5;
8479 n_floor = (
G4int)(r_floor);
8480 if (n_floor%2 == 0) {
8482 r_help = r_middle + (r_rest - 0.5) * (1.0 - r_even_odd);
8486 r_help = r_middle + (r_rest - 0.5) * (1.0 + r_even_odd);
8488 i_out = (
G4int)(r_help);
8504 G4double xcom = 0.0, xvs = 0.0, xe = 0.0;
8508 alpha = ( std::sqrt(5.0/(4.0*pi)) ) * beta;
8510 xcom = 1.0 - 1.7826 * ((a - 2.0*
z)/a)*((a - 2.0*
z)/a);
8512 xvs = - xcom * ( 15.4941 * a -
8513 17.9439 * std::pow(a,2.0/3.0) * (1.0+0.4*alpha*
alpha) );
8515 xe = z*z * (0.7053/(std::pow(a,1.0/3.0)) * (1.0-0.2*alpha*alpha) - 1.1529/
a);
8542 dtot = r0 * ( std::pow((z1+n1),1.0/3.0) * (1.0+0.6666667*beta1)
8543 + std::pow((z2+n2),1.0/3.0) * (1.0+0.6666667*beta2) ) +
d;
8544 fecoul = z1 * z2 * 1.44 / dtot;
8556 R_wash = std::exp(-E * Freduction * gamma);
8558 R_wash = std::exp(- Ecrit * Freduction * gamma -(E-Ecrit) * gamma);
8616 G4double z = 0.0,
n = 0.0,
a = 0.0, av = 0.0, as = 0.0;
8617 G4double a0 = 0.0,
c1 = 0.0, c4 = 0.0, b1 = 0.0, b3 = 0.0;
8619 G4double r0 = 0.0, kf = 0.0, ks = 0.0;
8620 G4double kv = 0.0, rp = 0.0, ay = 0.0, aden = 0.0, x0 = 0.0, y0 = 0.0;
8621 G4double esq = 0.0, ael = 0.0, i = 0.0;
8672 c1 = 3.0/5.0*esq/r0;
8673 c4 = 5.0/4.0*std::pow((3.0/(2.0*pi)),(2.0/3.0)) *
c1;
8674 kf = std::pow((9.0*pi*z/(4.0*
a)),(1.0/3.0))/r0;
8676 ff = -1.0/8.0*rp*rp*esq/std::pow(r0,3) * (145.0/48.0 - 327.0/2880.0*std::pow(kf,2) * std::pow(rp,2) + 1527.0/1209600.0*std::pow(kf,4) * std::pow(rp,4));
8680 x0 = r0 * std::pow(
a,(1.0/3.0)) / ay;
8681 y0 = r0 * std::pow(
a,(1.0/3.0)) / aden;
8683 b1 = 1.0 - 3.0/(std::pow(x0,2)) + (1.0 + x0) * (2.0 + 3.0/x0 + 3.0/std::pow(x0,2)) * std::exp(-2.0*x0);
8685 b3 = 1.0 - 5.0/std::pow(y0,2) * (1.0 - 15.0/(8.0*y0) + 21.0/(8.0 * std::pow(y0,3))
8686 - 3.0/4.0 * (1.0 + 9.0/(2.0*y0) + 7.0/std::pow(y0,2)
8687 + 7.0/(2.0 * std::pow(y0,3))) * std::exp(-2.0*y0));
8691 efl = -1.0 * av*(1.0 - kv*i*i)*
a + as*(1.0 - ks*i*i)*b1 * std::pow(
a,(2.0/3.0)) + a0
8692 +
c1*z*z*b3/std::pow(
a,(1.0/3.0)) - c4*std::pow(z,(4.0/3.0))/std::pow(
a,(1.e0/3.e0))
8693 + ff*std::pow(z,2)/
a -ca*(
n-
z) - ael * std::pow(z,(2.39e0));
8699 return eflmacResult;
8704 void G4Abla::unstable_nuclei(
G4int AFP,
G4int ZFP,
G4int *AFPNEW,
G4int *ZFPNEW,
G4int &IOUNSTABLE,
G4double VX,
G4double VY,
G4double VZ,
G4double *VP1X,
G4double *VP1Y,
G4double *VP1Z,
G4double BU_TAB_TEMP[200][5],
G4int *ILOOP){
8706 G4int INMIN,INMAX,NDIF=0,IMEM;
8707 G4int NEVA=0,PEVA=0;
8715 for(
G4int i=0;i<200;i++){
8716 BU_TAB_TEMP[i][0] = 0.0;
8717 BU_TAB_TEMP[i][1] = 0.0;
8718 BU_TAB_TEMP[i][2] = 0.0;
8719 BU_TAB_TEMP[i][3] = 0.0;
8720 BU_TAB_TEMP[i][4] = 0.0;
8726 if(AFP==0 && ZFP==0){
8730 if((AFP==1 && ZFP==0) ||
8731 (AFP==1 && ZFP==1) ||
8732 (AFP==2 && ZFP==1) ||
8733 (AFP==3 && ZFP==1) ||
8734 (AFP==3 && ZFP==2) ||
8735 (AFP==4 && ZFP==2) ||
8736 (AFP==6 && ZFP==2) ||
8745 if ((AFP-ZFP)==0 && ZFP>1){
8746 for(
G4int I = 0;I<=AFP-2;I++){
8747 unstable_tke(
double(AFP-I),
double(AFP-I),
double(AFP-I-1),
double(AFP-I-1),VX,VY,VZ,
8748 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
8749 BU_TAB_TEMP[*ILOOP][0] = 1.0;
8750 BU_TAB_TEMP[*ILOOP][1] = 1.0;
8751 BU_TAB_TEMP[*ILOOP][2] = VP2X;
8752 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
8753 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
8754 *ILOOP = *ILOOP + 1;
8759 PEVA = PEVA + ZFP - 1;
8771 for(
G4int I = 1;I<=10; I++){
8783 for(
G4int I = 0;I< IMEM;I++){
8789 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
8790 BU_TAB_TEMP[I+1+*ILOOP][0] = 1.0;
8791 BU_TAB_TEMP[I+1+*ILOOP][1] = 1.0;
8792 BU_TAB_TEMP[I+1+*ILOOP][2] = VP2X;
8793 BU_TAB_TEMP[I+1+*ILOOP][3] = VP2Y;
8794 BU_TAB_TEMP[I+1+*ILOOP][4] = VP2Z;
8799 *ILOOP = *ILOOP + IMEM;
8804 NEVA = NDIF - INMAX;
8807 for(
G4int I = 0;I<NEVA;I++){
8813 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
8815 BU_TAB_TEMP[*ILOOP][0] = 0.0;
8816 BU_TAB_TEMP[*ILOOP][1] = 1.0;
8817 BU_TAB_TEMP[*ILOOP][2] = VP2X;
8818 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
8819 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
8820 *ILOOP = *ILOOP + 1;
8827 if ((AFP>=2) && (ZFP==0)){
8828 for(
G4int I = 0;I<= AFP-2;I++){
8832 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
8834 BU_TAB_TEMP[*ILOOP][0] = 0.0;
8835 BU_TAB_TEMP[*ILOOP][1] = 1.0;
8836 BU_TAB_TEMP[*ILOOP][2] = VP2X;
8837 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
8838 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
8839 *ILOOP = *ILOOP + 1;
8845 NEVA = NEVA + (AFP - 1);
8851 std::cout <<
"WARNING - BU UNSTABLE: AF < ZF" << std::endl;
8856 if ((AFP>=4) && (ZFP==1)){
8858 for(
G4int I = 0; I<AFP-3;I++){
8860 double(AFP-I-1),
double(ZFP),
8862 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
8864 BU_TAB_TEMP[*ILOOP][0] = 0.0;
8865 BU_TAB_TEMP[*ILOOP][1] = 1.0;
8866 BU_TAB_TEMP[*ILOOP][2] = VP2X;
8867 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
8868 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
8869 *ILOOP = *ILOOP + 1;
8875 NEVA = NEVA + (AFP - 3);
8881 if ((AFP==4) && (ZFP==3)){
8889 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
8891 BU_TAB_TEMP[*ILOOP][0] = 1.0;
8892 BU_TAB_TEMP[*ILOOP][1] = 1.0;
8893 BU_TAB_TEMP[*ILOOP][2] = VP2X;
8894 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
8895 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
8896 *ILOOP = *ILOOP + 1;
8898 if ((AFP==5) && (ZFP==2)){
8906 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
8907 BU_TAB_TEMP[*ILOOP][0] = 0.0;
8908 BU_TAB_TEMP[*ILOOP][1] = 1.0;
8909 BU_TAB_TEMP[*ILOOP][2] = VP2X;
8910 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
8911 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
8912 *ILOOP = *ILOOP + 1;
8915 if ((AFP==5) && (ZFP==3)){
8923 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
8924 BU_TAB_TEMP[*ILOOP][0] = 1.0;
8925 BU_TAB_TEMP[*ILOOP][1] = 1.0;
8926 BU_TAB_TEMP[*ILOOP][2] = VP2X;
8927 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
8928 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
8929 *ILOOP = *ILOOP + 1;
8932 if ((AFP==6) && (ZFP==4)){
8941 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
8942 BU_TAB_TEMP[*ILOOP][0] = 1.0;
8943 BU_TAB_TEMP[*ILOOP][1] = 1.0;
8944 BU_TAB_TEMP[*ILOOP][2] = VP2X;
8945 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
8946 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
8947 *ILOOP = *ILOOP + 1;
8955 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
8956 BU_TAB_TEMP[*ILOOP][0] = 1.0;
8957 BU_TAB_TEMP[*ILOOP][1] = 1.0;
8958 BU_TAB_TEMP[*ILOOP][2] = VP2X;
8959 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
8960 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
8961 *ILOOP = *ILOOP + 1;
8963 if ((AFP==7)&&(ZFP==2)){
8971 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
8972 BU_TAB_TEMP[*ILOOP][0] = 0.0;
8973 BU_TAB_TEMP[*ILOOP][1] = 1.0;
8974 BU_TAB_TEMP[*ILOOP][2] = VP2X;
8975 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
8976 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
8977 *ILOOP = *ILOOP + 1;
8980 if ((AFP==7) && (ZFP==5)){
8982 for(
G4int I = 0; I<= AFP-5;I++){
8984 double(AFP-I-1),
double(ZFP-I-1),
8986 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
8987 BU_TAB_TEMP[*ILOOP][0] = 1.0;
8988 BU_TAB_TEMP[*ILOOP][1] = 1.0;
8989 BU_TAB_TEMP[*ILOOP][2] = VP2X;
8990 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
8991 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
8992 *ILOOP = *ILOOP + 1;
9003 if ((AFP==8) && (ZFP==4)){
9010 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9011 BU_TAB_TEMP[*ILOOP][0] = 2.0;
9012 BU_TAB_TEMP[*ILOOP][1] = 4.0;
9013 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9014 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9015 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9016 *ILOOP = *ILOOP + 1;
9018 if ((AFP==8) && (ZFP==6)){
9027 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9028 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9029 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9030 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9031 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9032 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9033 *ILOOP = *ILOOP + 1;
9040 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9041 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9042 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9043 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9044 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9045 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9046 *ILOOP = *ILOOP + 1;
9052 if((AFP==9) && (ZFP==2)){
9061 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9062 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9063 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9064 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9065 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9066 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9067 *ILOOP = *ILOOP + 1;
9073 if((AFP==9) && (ZFP==5)){
9081 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9082 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9083 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9084 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9085 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9086 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9087 *ILOOP = *ILOOP + 1;
9094 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9095 BU_TAB_TEMP[*ILOOP][0] = 2.0;
9096 BU_TAB_TEMP[*ILOOP][1] = 4.0;
9097 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9098 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9099 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9100 *ILOOP = *ILOOP + 1;
9106 if((AFP==10) && (ZFP==2)){
9115 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9116 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9117 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9118 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9119 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9120 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9121 *ILOOP = *ILOOP + 1;
9129 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9130 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9131 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9132 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9133 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9134 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9135 *ILOOP = *ILOOP + 1;
9140 if ((AFP==10) && (ZFP==3)){
9148 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9149 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9150 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9151 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9152 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9153 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9154 *ILOOP = *ILOOP + 1;
9159 if ((AFP==10) && (ZFP==7)){
9167 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9168 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9169 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9170 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9171 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9172 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9173 *ILOOP = *ILOOP + 1;
9179 if((AFP==11) && (ZFP==7)){
9187 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9188 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9189 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9190 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9191 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9192 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9193 *ILOOP = *ILOOP + 1;
9198 if ((AFP==12) && (ZFP==8)){
9207 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9208 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9209 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9210 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9211 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9212 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9213 *ILOOP = *ILOOP + 1;
9220 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9221 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9222 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9223 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9224 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9225 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9226 *ILOOP = *ILOOP + 1;
9231 if ((AFP==15) && (ZFP==9)){
9239 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9240 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9241 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9242 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9243 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9244 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9245 *ILOOP = *ILOOP + 1;
9251 if ((AFP==16) && (ZFP==9)){
9259 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9260 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9261 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9262 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9263 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9264 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9265 *ILOOP = *ILOOP + 1;
9271 if ((AFP==16) && (ZFP==10)){
9279 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9280 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9281 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9282 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9283 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9284 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9285 *ILOOP = *ILOOP + 1;
9292 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9293 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9294 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9295 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9296 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9297 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9298 *ILOOP = *ILOOP + 1;
9303 if((AFP==18) && (ZFP==11)){
9311 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9312 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9313 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9314 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9315 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9316 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9317 *ILOOP = *ILOOP + 1;
9322 if((AFP==19) && (ZFP==11)){
9330 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9331 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9332 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9333 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9334 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9335 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9336 *ILOOP = *ILOOP + 1;
9341 if (ZFP>=4 && (AFP-ZFP)==1){
9346 for(
G4int I = 0;I< NEVA;I++){
9350 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9351 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9352 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9353 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9354 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9355 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9356 *ILOOP = *ILOOP + 1;
9361 for(
G4int I = 0;I<PEVA;I++){
9365 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9366 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9367 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9368 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9369 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9370 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9371 *ILOOP = *ILOOP + 1;
9389 void G4Abla::unstable_tke(
G4double ain,
G4double zin,
G4double anew,
G4double znew,
G4double vxin,
G4double vyin,
G4double vzin,
G4double *v1x,
G4double *v1y,
G4double *v1z,
G4double *v2x,
G4double *v2y,
G4double *v2z){
9392 G4double PX1,PX2,PY1,PY2,PZ1,PZ2,PTOT;
9393 G4double RNDT,CTET1,STET1,RNDP,PHI1,ETOT_P1,ETOT_P2;
9395 G4double vxout=0.,vyout=0.,vzout=0.;
9396 G4int iain,izin,ianew,iznew,inin,innew;
9406 innew = ianew - iznew;
9409 mglms(ain,zin,3,&MASS);
9410 mglms(anew,znew,3,&MASS1);
9411 mglms(ain-anew,zin-znew,3,&MASS2);
9412 ekin_tot = MASS-MASS1-MASS2;
9416 if(izin>12)std::cout <<
"*** ZIN > 12 ***" << izin << std::endl;
9419 if( ekin_tot<0.00 ){
9430 EKIN_P1 = ekin_tot*(ain-anew)/ ain;
9431 ETOT_P1 = EKIN_P1 + anew * AMU;
9432 PTOT = anew*AMU*std::sqrt((EKIN_P1/(anew*AMU)+1.0)*(EKIN_P1/(anew*AMU)+1.0)-1.0);
9435 CTET1 = 2.0*RNDT-1.0;
9436 STET1 = std::sqrt(1.0-CTET1*CTET1);
9438 PHI1 = RNDP*2.0*3.141592654;
9439 PX1 = PTOT * STET1*std::cos(PHI1);
9440 PY1 = PTOT * STET1*std::sin(PHI1);
9442 *v1x = C * PX1 / ETOT_P1;
9443 *v1y = C * PY1 / ETOT_P1;
9444 *v1z = C * PZ1 / ETOT_P1;
9445 lorentz_boost(vxin,vyin,vzin,*v1x,*v1y,*v1z,&vxout,&vyout,&vzout);
9453 ETOT_P2 = (ekin_tot - EKIN_P1) + (ain-anew) * AMU;
9454 *v2x = C * PX2 / ETOT_P2;
9455 *v2y = C * PY2 / ETOT_P2;
9456 *v2z = C * PZ2 / ETOT_P2;
9457 lorentz_boost(vxin,vyin,vzin,*v2x,*v2y,*v2z,&vxout,&vyout,&vzout);
9475 G4double GAMMA,VR,
C,CC,DENO,VXNOM,VYNOM,VZNOM;
9486 VR = std::sqrt(VXR*VXR + VYR*VYR + VZR*VZR);
9493 GAMMA = 1.0/std::sqrt(1.0 - VR*VR/CC);
9494 DENO = 1.0 - VXR*VXIN/CC - VYR*VYIN/CC - VZR*VZIN/CC;
9497 VXNOM = -GAMMA*VXR + (1.0+(GAMMA-1.0)*VXR*VXR/(VR*VR))*VXIN + (GAMMA-1.0)*VXR*VYR/(VR*VR)*VYIN + (GAMMA-1.0)*VXR*VZR/(VR*VR)*VZIN;
9499 *VXOUT = VXNOM / (GAMMA * DENO);
9502 VYNOM = -GAMMA*VYR + (1.0+(GAMMA-1.0)*VYR*VYR/(VR*VR))*VYIN + (GAMMA-1.0)*VXR*VYR/(VR*VR)*VXIN + (GAMMA-1.0)*VYR*VZR/(VR*VR)*VZIN;
9504 *VYOUT = VYNOM / (GAMMA * DENO);
9507 VZNOM = -GAMMA*VZR + (1.0+(GAMMA-1.0)*VZR*VZR/(VR*VR))*VZIN + (GAMMA-1.0)*VXR*VZR/(VR*VR)*VXIN + (GAMMA-1.0)*VYR*VZR/(VR*VR)*VYIN;
9509 *VZOUT = VZNOM / (GAMMA * DENO);
9521 G4double EFF1=0.,EFF2=0.,VFF1=0.,VFF2=0.,
9522 AF1=0.,ZF1=0.,AF2=0.,ZF2=0.,
9523 AFF1=0.,ZFF1=0.,AFF2=0.,ZFF2=0.,
9524 vz1_eva=0., vx1_eva=0.,vy1_eva=0.,
9525 vz2_eva=0., vx2_eva=0.,vy2_eva=0.,
9526 vx_eva_sc=0.,vy_eva_sc=0.,vz_eva_sc=0.,
9527 VXOUT=0.,VYOUT=0.,VZOUT=0.,
9528 VX2OUT=0.,VY2OUT=0.,VZ2OUT=0.;
9529 G4int IEV_TAB_FIS=0,IEV_TAB_TEMP=0;
9530 G4double EV_TEMP1[200][5], EV_TEMP2[200][5],mtota;
9531 G4int inttype = 0,inum=0;
9535 for(
G4int I1=0;I1<200;I1++)
9536 for(
G4int I2=0;I2<5;I2++){
9537 EV_TEMP[I1][I2] = 0.0;
9538 EV_TEMP1[I1][I2] = 0.0;
9539 EV_TEMP2[I1][I2] = 0.0;
9542 G4double et = EE - JPRF * JPRF * 197. * 197./(2.*0.4*931.*std::pow(AF,5.0/3.0)*1.16*1.16);
9544 fissionDistri(AF,ZF,et,AF1,ZF1,EFF1,VFF1,AF2,ZF2,EFF2,VFF2,
9545 vx_eva_sc,vy_eva_sc,vz_eva_sc);
9559 G4double VPERP1 = std::sqrt(VFF1*VFF1 - VZ1_FISSION*VZ1_FISSION);
9561 G4double VX1_FISSION = VPERP1 * std::sin(ALPHA1);
9562 G4double VY1_FISSION = VPERP1 * std::cos(ALPHA1);
9563 G4double VX2_FISSION = - VX1_FISSION / VFF1 * VFF2;
9564 G4double VY2_FISSION = - VY1_FISSION / VFF1 * VFF2;
9565 G4double VZ2_FISSION = - VZ1_FISSION / VFF1 * VFF2;
9568 if( (ZF1<=0.0) || (AF1<=0.0) || (AF1<ZF1) ){
9569 std::cout <<
"F1 unphysical: "<<ZF<<
" "<<AF<<
" "<<EE<<
" "<<ZF1<<
" "<<AF1 << std::endl;
9575 G4int FF11=0, FIMF11=0;
9576 G4double ZIMFF1=0., AIMFF1=0.,TKEIMF1=0.,JPRFOUT=0.;
9578 evapora(ZF1,AF1,&EFF1,0., &ZFF1, &AFF1, &mtota, &vz1_eva, &vx1_eva,&vy1_eva, &FF11, &FIMF11, &ZIMFF1, &AIMFF1,&TKEIMF1, &JPRFOUT, &inttype, &inum,EV_TEMP1,&IEV_TAB_TEMP);
9580 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
9581 EV_TEMP[IJ+IEV_TAB_FIS][0] = EV_TEMP1[IJ][0];
9582 EV_TEMP[IJ+IEV_TAB_FIS][1] = EV_TEMP1[IJ][1];
9589 EV_TEMP1[IJ][2],EV_TEMP1[IJ][3],EV_TEMP1[IJ][4],
9590 &VXOUT,&VYOUT,&VZOUT);
9593 &VX2OUT,&VY2OUT,&VZ2OUT);
9594 EV_TEMP[IJ+IEV_TAB_FIS][2] = VX2OUT;
9595 EV_TEMP[IJ+IEV_TAB_FIS][3] = VY2OUT;
9596 EV_TEMP[IJ+IEV_TAB_FIS][4] = VZ2OUT;
9599 IEV_TAB_FIS = IEV_TAB_FIS + IEV_TAB_TEMP;
9604 if( (ZF2<=0.0) || (AF2<=0.0) || (AF2<ZF2) ){
9605 std::cout <<
"F2 unphysical: "<<ZF<<
" "<<AF<<
" "<<EE<<
" "<<ZF2<<
" "<<AF2 << std::endl;
9611 G4int FF22=0, FIMF22=0;
9612 G4double ZIMFF2=0., AIMFF2=0.,TKEIMF2=0.,JPRFOUT=0.;
9614 evapora(ZF2,AF2,&EFF2,0., &ZFF2, &AFF2, &mtota, &vz2_eva, &vx2_eva,&vy2_eva, &FF22, &FIMF22, &ZIMFF2, &AIMFF2,&TKEIMF2, &JPRFOUT, &inttype, &inum,EV_TEMP2,&IEV_TAB_TEMP);
9616 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
9617 EV_TEMP[IJ+IEV_TAB_FIS][0] = EV_TEMP2[IJ][0];
9618 EV_TEMP[IJ+IEV_TAB_FIS][1] = EV_TEMP2[IJ][1];
9625 EV_TEMP2[IJ][2],EV_TEMP2[IJ][3],EV_TEMP2[IJ][4],
9626 &VXOUT,&VYOUT,&VZOUT);
9629 &VX2OUT,&VY2OUT,&VZ2OUT);
9630 EV_TEMP[IJ+IEV_TAB_FIS][2] = VX2OUT;
9631 EV_TEMP[IJ+IEV_TAB_FIS][3] = VY2OUT;
9632 EV_TEMP[IJ+IEV_TAB_FIS][4] = VZ2OUT;
9635 IEV_TAB_FIS = IEV_TAB_FIS + IEV_TAB_TEMP;
9648 VX1_FISSION,VY1_FISSION,VZ1_FISSION,
9649 &VXOUT,&VYOUT,&VZOUT);
9650 VX1_FISSION = VXOUT;
9651 VY1_FISSION = VYOUT;
9652 VZ1_FISSION = VZOUT;
9654 VX2_FISSION,VY2_FISSION,VZ2_FISSION,
9655 &VXOUT,&VYOUT,&VZOUT);
9656 VX2_FISSION = VXOUT;
9657 VY2_FISSION = VYOUT;
9658 VZ2_FISSION = VZOUT;
9662 (*VX1_FISSION_par) = VX1_FISSION;
9663 (*VY1_FISSION_par) = VY1_FISSION;
9664 (*VZ1_FISSION_par) = VZ1_FISSION;
9665 (*VX_EVA_SC_par)=vx_eva_sc;
9666 (*VY_EVA_SC_par)=vy_eva_sc;
9667 (*VZ_EVA_SC_par)=vz_eva_sc;
9670 (*VX2_FISSION_par) = VX2_FISSION;
9671 (*VY2_FISSION_par) = VY2_FISSION;
9672 (*VZ2_FISSION_par) = VZ2_FISSION;
9673 (*IEV_TAB_FIS_par) = IEV_TAB_FIS;
9680 G4double V_over_V0,R0,RALL,RHAZ,
R,TKE,Ekin,V,VPERP,ALPHA1;
9692 RALL = R0 * std::pow(V_over_V0,1.0/3.0) * std::pow(AAL,1.0/3.0);
9694 R = std::pow(RHAZ,1.0/3.0) * RALL;
9695 TKE = 1.44 * Z * ZALL * R*R * (1.0 - A/AAL)*(1.0 - A/AAL) / std::pow(RALL,3.0);
9697 Ekin = TKE * (AAL -
A) / AAL;
9699 V = std::sqrt(Ekin/A) * 1.3887;
9701 VPERP = std::sqrt(V*V - (*VZ)*(*VZ));
9703 *VX = VPERP * std::sin(ALPHA1);
9704 *VY = VPERP * std::cos(ALPHA1);
9730 ix =
G4int(
y * 100 + 43543000);
9731 if(
mod(ix,2) == 0) {
9745 l_plus = lambda + 1.;
9749 y=std::pow(
G4AblaRandom::flat()*(std::pow(rxmax,l_plus)-std::pow(rxmin,l_plus))+ std::pow(rxmin,l_plus),1.0/l_plus);
9761 G4double Efermi = 5.0 * SIGMA_0 * SIGMA_0 / (2.0 * 931.4940);
9767 SIGMA_0 = SIGMA_0 * std::pow(V0_over_VBU,1.0/3.0);
9772 GOLDHA_BU = SIGMA_0 * std::sqrt((APRF*(AABRA-APRF))/(AABRA-1.0));
9773 GOLDHA = GOLDHA_BU*std::sqrt(1.0 +
9781 (AABRA - APRF) / AABRA);
9787 GOLDHA = SIGMA_0 * std::sqrt((APRF*(AABRA-APRF))/(AABRA-1.0));
9795 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING PX IN Rn07.FOR. A VALUE WILL BE FORCED." << std::endl;
9796 *PX = (AABRA-1.0)*931.4940;
9798 if(std::abs(*PX)>= AABRA*931.494){
9807 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING PY IN Rn07.FOR. A VALUE WILL BE FORCED." << std::endl;
9808 *PY = (AABRA-1.0)*931.4940;
9810 if(std::abs(*PY)>= AABRA*931.494){
9819 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING PZ IN Rn07.FOR. A VALUE WILL BE FORCED." << std::endl;
9820 *PZ = (AABRA-1.0)*931.4940;
9822 if(std::abs(*PZ)>= AABRA*931.494){
9839 v1 = 2.0*
haz(k) - 1.0;
9840 v2 = 2.0*
haz(k) - 1.0;
9841 r = std::pow(v1,2) + std::pow(v2,2);
9844 fac = std::sqrt(-2.*std::log(r)/r);
9846 fgausshaz = v2*fac*sig+xmoy;
9850 fgausshaz=gset*sig+xmoy;
void AMOMENT(G4double AABRA, G4double APRF, G4int IMULTIFR, G4double *PX, G4double *PY, G4double *PZ)
G4double fmaxhaz_old(G4double T)
void FillData(G4int IMULTBU, G4int IEV_TAB)
G4double gausshaz(int k, double xmoy, double sig)
std::vector< ExP01TrackerHit * > a
void part_fiss(G4double BET, G4double GP, G4double GF, G4double Y, G4double TAUF, G4double TS1, G4double TSUM, G4int *CHOICE, G4double ZF, G4double AF, G4double FT, G4double *T_LAPSE, G4double *GF_LOC)
G4double tau(G4double bet, G4double homega, G4double ef, G4double t)
G4double getMexp(G4int A, G4int Z)
void unbound(G4double SN, G4double SP, G4double SD, G4double ST, G4double SHE, G4double SA, G4double BP, G4double BD, G4double BT, G4double BHE, G4double BA, G4double *PROBF, G4double *PROBN, G4double *PROBP, G4double *PROBD, G4double *PROBT, G4double *PROBHE, G4double *PROBA, G4double *PROBIMF, G4double *PROBG, G4double *ECN, G4double *ECP, G4double *ECD, G4double *ECT, G4double *ECHE, G4double *ECA)
G4double tunnelling(G4double A, G4double ZPRF, G4double Y, G4double EE, G4double EF, G4double TEMP, G4double DENSG, G4double DENSF, G4double ENH_FACT)
void qrot(G4double z, G4double a, G4double bet, G4double sig, G4double u, G4double *qr)
Float_t y1[n_points_granero]
G4double massexp[MASSIZEROWS][MASSIZECOLS]
void mglw(G4double a, G4double z, G4double *el)
void mglms(G4double a, G4double z, G4int refopt4, G4double *el)
Float_t x1[n_points_granero]
void SetParametersG4(G4int z, G4int a)
void appariem(G4double a, G4double z, G4double *del)
void fissionDistri(G4double &a, G4double &z, G4double &e, G4double &a1, G4double &z1, G4double &e1, G4double &v1, G4double &a2, G4double &z2, G4double &e2, G4double &v2, G4double &vx_eva_sc, G4double &vy_eva_sc, G4double &vz_eva_sc)
void bsbkbc(G4double A, G4double Z, G4double *BS, G4double *BK, G4double *BC)
G4double umass(G4double z, G4double n, G4double beta)
void DeexcitationAblaxx(G4int nucleusA, G4int nucleusZ, G4double excitationEnergy, G4double angularMomentum, G4double momX, G4double momY, G4double momZ, G4int eventnumber)
void guet(G4double *x_par, G4double *z_par, G4double *find_par)
G4double pen(G4double A, G4double ap, G4double omega, G4double T)
G4double getEcnz(G4int A, G4int Z)
G4int idnint(G4double value)
void even_odd(G4double r_origin, G4double r_even_odd, G4int &i_out)
G4int nint(G4double number)
void unstable_tke(G4double AIN, G4double ZIN, G4double ANEW, G4double ZNEW, G4double VXIN, G4double VYIN, G4double VZIN, G4double *V1X, G4double *V1Y, G4double *V1Z, G4double *V2X, G4double *V2Y, G4double *V2Z)
G4double cram(G4double bet, G4double homega)
void fission(G4double AF, G4double ZF, G4double EE, G4double JPRF, G4double *VX1_FISSION, G4double *VY1_FISSION, G4double *VZ1_FISSION, G4double *VX2_FISSION, G4double *VY2_FISSION, G4double *VZ2_FISSION, G4int *ZFP1, G4int *AFP1, G4int *ZFP2, G4int *AFP2, G4int *imode, G4double *VX_EVA_SC, G4double *VY_EVA_SC, G4double *VZ_EVA_SC, G4double EV_TEMP[200][5], G4int *IEV_TAB_FIS)
G4double ecoul(G4double z1, G4double n1, G4double beta1, G4double z2, G4double n2, G4double beta2, G4double d)
G4double func_trans(G4double TIME, G4double ZF, G4double AF, G4double BET, G4double Y, G4double FT, G4double T_0)
Float_t y2[n_points_geant4]
void lpoly(G4double x, G4int n, G4double pl[])
G4double beta2[ECLDROWSbeta][ECLDCOLSbeta]
G4double dint(G4double a)
static constexpr double m
static constexpr double pc
G4double ecgnz[ECLDROWS][ECLDCOLS]
G4int getMexpID(G4int A, G4int Z)
void evap_postsaddle(G4double A, G4double Z, G4double E_scission_pre, G4double *E_scission_post, G4double *A_scission, G4double *Z_scission, G4double &vx_eva, G4double &vy_eva, G4double &vz_eva)
G4double eflmac(G4int ia, G4int iz, G4int flag, G4int optshp)
G4double enerj[VARNTPSIZE]
G4double pzlab[VARNTPSIZE]
G4double getAlpha(G4int A, G4int Z)
const XML_Char int const XML_Char * value
G4int min(G4int a, G4int b)
G4double pylab[VARNTPSIZE]
G4double frldm(double z, double n, double beta)
G4int mod(G4int a, G4int b)
G4double getBeta2(G4int A, G4int Z)
G4double gammp(G4double a, G4double x)
static const G4double alpha
G4int itypcasc[VARNTPSIZE]
void fission_width(G4double ZPRF, G4double A, G4double EE, G4double BS, G4double BK, G4double EF, G4double Y, G4double *GF, G4double *TEMP, G4double JPR, G4int IEROT, G4int FF_ALLOWED, G4int OPTCOL, G4int OPTSHP, G4double DENSG)
double A(double temperature)
void barrs(G4int Z1, G4int A1, G4int Z2, G4int A2, G4double *sBARR, G4double *sOMEGA)
G4int ISIGN(G4int a, G4int b)
void fomega_sp(G4double AF, G4double Y, G4double *MFCD, G4double *sOMEGA, G4double *sHOMEGA)
G4double pxlab[VARNTPSIZE]
void gcf(G4double *gammcf, G4double a, G4double x, G4double gln)
G4double alpha[ECLDROWS][ECLDCOLS]
G4double dmin1(G4double a, G4double b, G4double c)
G4double getRms(G4int A, G4int Z)
void setVerboseLevel(G4int level)
void densniv(G4double a, G4double z, G4double ee, G4double ef, G4double *dens, G4double bshell, G4double bs, G4double bk, G4double *temp, G4int optshp, G4int optcol, G4double defbet, G4double *ecor, G4double jprf, G4int ifis, G4double *qr)
G4double utilabs(G4double a)
G4double rms[ECLDROWS][ECLDCOLS]
G4double DSIGN(G4double a, G4double b)
void parite(G4double n, G4double *par)
G4double eflmac_profi(double a, double z)
void imf(G4double ACN, G4double ZCN, G4double TEMP, G4double EE, G4double *ZIMF, G4double *AIMF, G4double *BIMF, G4double *SBIMF, G4double *TIMF, G4double JPRF)
G4Abla(G4Volant *aVolant, G4VarNtp *aVarntp)
G4double spdef(G4int a, G4int z, G4int optxfis)
G4double vgsld[ECLDROWS][ECLDCOLS]
G4double getBeta4(G4int A, G4int Z)
G4double dm[PACESIZEROWS][PACESIZECOLS]
void lorentz_boost(G4double VXRIN, G4double VYRIN, G4double VZRIN, G4double VXIN, G4double VYIN, G4double VZIN, G4double *VXOUT, G4double *VYOUT, G4double *VZOUT)
G4int max(G4int a, G4int b)
void isostab_lim(G4int z, G4int *nmin, G4int *nmax)
G4int IPOWERLIMHAZ(G4double lambda, G4int xmin, G4int xmax)
static const G4double fac
G4double fvmaxhaz_neut(G4double x)
G4double getPace2(G4int A, G4int Z)
G4double EV_TAB_SSC[200][5]
void fomega_gs(G4double AF, G4double ZF, G4double *K1, G4double *sOMEGA, G4double *sHOMEGA)
void direct(G4double zprf, G4double a, G4double ee, G4double jprf, G4double *probp_par, G4double *probd_par, G4double *probt_par, G4double *probn_par, G4double *probhe_par, G4double *proba_par, G4double *probg_par, G4double *probimf_par, G4double *probf_par, G4double *ptotl_par, G4double *sn_par, G4double *sbp_par, G4double *sbd_par, G4double *sbt_par, G4double *sbhe_par, G4double *sba_par, G4double *ecn_par, G4double *ecp_par, G4double *ecd_par, G4double *ect_par, G4double *eche_par, G4double *eca_par, G4double *ecg_par, G4double *bp_par, G4double *bd_par, G4double *bt_par, G4double *bhe_par, G4double *ba_par, G4double *sp, G4double *sd, G4double *st, G4double *she, G4double *sa, G4double *ef, G4double *ts1, G4int inttype, G4int inum, G4int itest, G4int *sortie, G4double *tcn, G4double *jprfn, G4double *jprfp, G4double *jprfd, G4double *jprft, G4double *jprfhe, G4double *jprfa, G4double *tsum)
G4double gammln(G4double xx)
G4double pace2(G4double a, G4double z)
void gser(G4double *gamser, G4double a, G4double x, G4double gln)
G4double fvmaxhaz(G4double T)
void evapora(G4double zprf, G4double aprf, G4double *ee_par, G4double jprf, G4double *zf_par, G4double *af_par, G4double *mtota_par, G4double *vleva_par, G4double *vxeva_par, G4double *vyeva_par, G4int *ff_par, G4int *fimf_par, G4double *fzimf, G4double *faimf, G4double *tkeimf_par, G4double *jprfout, G4int *inttype_par, G4int *inum_par, G4double EV_TEMP[200][5], G4int *iev_tab_temp_par)
static constexpr double pi
static constexpr double amu
void unstable_nuclei(G4int AFP, G4int ZFP, G4int *AFPNEW, G4int *ZFPNEW, G4int &IOUNSTABLE, G4double VX, G4double VY, G4double VZ, G4double *VP1X, G4double *VP1Y, G4double *VP1Z, G4double BU_TAB_TEMP[200][5], G4int *ILOOP)
static const G4double eps
G4double efa[FBCOLS][FBROWS]
G4double fmaxhaz(G4double T)
G4double ecfnz[ECLDROWS][ECLDCOLS]
void lorb(G4double AMOTHER, G4double ADAUGHTER, G4double LMOTHER, G4double EEFINAL, G4double *LORBITAL, G4double *SIGMA_LORBITAL)
G4double getVgsld(G4int A, G4int Z)
void tke_bu(G4double Z, G4double A, G4double ZALL, G4double AAL, G4double *VX, G4double *VY, G4double *VZ)
double B(double temperature)
G4double bind[MASSIZEROWS][MASSIZECOLS]
G4double bipol(int iflag, G4double y)
G4double fissility(int a, int z, int optxfis)
void barfit(G4int iz, G4int ia, G4int il, G4double *sbfis, G4double *segs, G4double *selmax)
G4int mexpiop[MASSIZEROWS][MASSIZECOLS]
G4double ecnz[EC2SUBROWS][EC2SUBCOLS]
G4double beta4[ECLDROWSbeta][ECLDCOLSbeta]
G4double width(G4double AMOTHER, G4double ZMOTHER, G4double APART, G4double ZPART, G4double TEMP, G4double B1, G4double SB1, G4double EXC)
G4double Uwash(double E, double Ecrit, double Freduction, double gamma)