34 #define ABLAXX_IN_GEANT4_MODE 1
43 #ifdef ABLAXX_IN_GEANT4_MODE
49 #ifndef ABLAXX_IN_GEANT4_MODE
98 DeexcitationAblaxx(nucleusA,nucleusZ,excitationEnergy,angularMomentum,momX,momY,momZ,eventnumber,0);
113 if(nucleusS>0)nucleusS=0;
123 G4int ZFP1 = 0, AFP1 = 0, AFPIMF = 0, ZFPIMF = 0, ZFP2 = 0, AFP2 = 0, SFP1 = 0, SFP2 = 0, SFPIMF = 0;
124 G4double vx_eva = 0.0, vy_eva = 0.0, vz_eva = 0.0;
125 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;
126 G4double zf = 0.0, af = 0.0, mtota = 0.0, tkeimf = 0.0, jprf0=0.;
127 G4int ff = 0,afpnew=0,zfpnew=0,aprfp=0,zprfp=0,IOUNSTABLE=0,ILOOP=0,IEV_TAB=0,IEV_TAB_TEMP=0;
128 G4int fimf = 0,INMIN=0,INMAX=0;
130 G4int inum = eventnumber;
160 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.;
162 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.;
164 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.;
166 G4int ABU_DIFF,ZBU_DIFF,NBU_DIFF;
167 G4int INEWLOOP = 0, ILOOPBU=0;
169 G4double BU_TAB_TEMP[200][6], BU_TAB_TEMP1[200][6];
170 G4double EV_TAB_TEMP[200][6],EV_TEMP[200][6];
171 G4int IMEM_BU[200], IMEM=0;
174 std::cout <<
"Error - Remnant with a mass number A below 1." << std::endl;
179 for(
G4int j=0;j<3;j++){
186 for(
G4int I1=0;I1<200;I1++){
187 for(
G4int I2 = 0;I2<12;I2++)
189 for(
G4int I2 = 0;I2<6;I2++){
190 BU_TAB_TEMP[I1][I2] = 0.0;
191 BU_TAB_TEMP1[I1][I2] = 0.0;
192 EV_TAB_TEMP[I1][I2] = 0.0;
195 EV_TEMP[I1][I2] = 0.0;
222 G4double pincl = std::sqrt(pxrem*pxrem + pyrem*pyrem + pzrem*pzrem);
224 G4double ETOT_incl = std::sqrt(pincl*pincl + (AAINCL * amu)*(AAINCL * amu));
225 G4double VX_incl = C * pxrem / ETOT_incl;
226 G4double VY_incl = C * pyrem / ETOT_incl;
227 G4double VZ_incl = C * pzrem / ETOT_incl;
265 a_tilda =
ald->
av*aprf +
ald->
as*std::pow(aprf,2.0/3.0) +
ald->
ak*std::pow(aprf,1.0/3.0);
267 T_init = std::sqrt(EINCL/a_tilda);
271 if(T_diff>0.1 && zprf>2. && (aprf-zprf)>0.){
276 for(
G4int i=0;i<5;i++){
283 A_diff =
dint(EE_diff / (8.0 * 5.0 / T_freeze_out));
285 if(A_diff>AAINCL) A_diff = AAINCL;
287 A_FINAL = AAINCL - A_diff;
289 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);
293 EE_diff = EINCL - E_FINAL;
305 Z_FINAL =
dint(zprf * A_FINAL / (aprf));
307 if(E_FINAL<0.0) E_FINAL = 0.0;
313 A_diff = AAINCL - aprf;
322 }
else if(A_diff>1.0){
330 a_tilda =
ald->
av*AAINCL +
ald->
as*std::pow(AAINCL,2.0/3.0) +
ald->
ak*std::pow(AAINCL,1.0/3.0);
334 ABU_SLOPE = (ASLOPE1-ASLOPE2)/4.0*(E_FINAL/AAINCL)+
335 ASLOPE1-(ASLOPE1-ASLOPE2)*7.0/4.0;
347 if(ABU_SLOPE > -1.01) ABU_SLOPE = -1.01;
350 Z_Breakup_sum = Z_FINAL;
354 for(
G4int i=0;i<100;i++){
361 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;
365 if(A_Breakup>AAINCL)
goto mult4326;
368 std::cout <<
"A_BREAKUP <= 0 " << std::endl;
372 A_ACC = A_ACC + A_Breakup;
376 Z_Breakup_Mean =
dint(A_Breakup * ZAINCL / AAINCL);
378 Z_Breakup_sum = Z_Breakup_sum + Z_Breakup_Mean;
381 G_SYMM = 34.2281 - 5.14037 * E_FINAL/AAINCL;
382 if(E_FINAL/AAINCL < 2.0) G_SYMM = 25.0;
383 if(E_FINAL/AAINCL > 4.0) G_SYMM = 15.0;
388 CZ = 2.0 * G_SYMM * 4.0 / A_Breakup;
392 Sigma_Z = std::sqrt(T_freeze_out/CZ);
400 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;
404 if(Z_Breakup<0.0 )
goto mult4333;
405 if((A_Breakup-Z_Breakup)<0.0)
goto mult4333;
406 if((A_Breakup-Z_Breakup)==0.0 && Z_Breakup!=1.0)
goto mult4333;
408 if(Z_Breakup>=ZAINCL){
411 std::cout <<
"Z_BREAKUP RESAMPLED MORE THAN 10 TIMES; EVENT WILL BE RESAMPLED AGAIN " << std::endl;
421 if(
idnint(A_Breakup-Z_Breakup)<INMIN ||
idnint(A_Breakup-Z_Breakup)>(INMAX+5)){
433 N_Breakup = A_Breakup - Z_Breakup;
436 ABU_SUM = ABU_SUM +
BU_TAB[i][1];
437 ZBU_SUM = ZBU_SUM + BU_TAB[i][0];
440 BU_TAB[I_Breakup][3] = 0.0;
441 I_Breakup = I_Breakup + 1;
442 IMULTBU = IMULTBU + 1;
459 ABU_DIFF =
idnint(ABU_SUM+aprf-AAINCL);
460 ZBU_DIFF =
idnint(ZBU_SUM+zprf-ZAINCL);
461 NBU_DIFF =
idnint((ABU_SUM-ZBU_SUM)+(aprf-zprf)-(AAINCL-ZAINCL));
464 std::cout <<
"WARNING - MORE THAN 200 BU " << IMULTBU << std::endl;
467 std::cout <<
"WARNING - LESS THAN 1 BU " << IMULTBU << std::endl;
471 for(
G4int i=0;i<IMULTBU;i++)
474 while(NBU_DIFF!=0 && ZBU_DIFF!=0){
483 std::cout <<
"WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED." << std::endl;
487 if(IMEM_BU[IEL]==1)
goto mult5555;
488 if(!(IEL<200))std::cout <<
"5555:" << IEL << RHAZ << IMULTBU << std::endl;
489 if(IEL<0)std::cout <<
"5555:"<< IEL << RHAZ << IMULTBU << std::endl;
492 }
else if(IEL>IMULTBU){
501 }
else if(IEL>IMULTBU){
508 if(ZTEMP<1.0 && N_Breakup<1.0){
522 }
else if(IEL>IMULTBU){
524 aprf =
dint(ZTEMP + N_Breakup);
526 NBU_DIFF = NBU_DIFF -
ISIGN(1,NBU_DIFF);
527 ZBU_DIFF = ZBU_DIFF -
ISIGN(1,ZBU_DIFF);
531 for(
G4int i=0;i<IMULTBU;i++)
534 if(NBU_DIFF != 0 && ZBU_DIFF == 0){
535 while(NBU_DIFF > 0 || NBU_DIFF < 0){
541 std::cout <<
"WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED." << std::endl;
545 if(IMEM_BU[IEL]==1)
goto mult5556;
547 if(IPROBA>IMULTBU+1 && NBU_DIFF>0){
548 std::cout <<
"###',IPROBA,IMULTBU,NBU_DIFF,ZBU_DIFF,T_freeze_out" << std::endl;
552 }
else{
if(IEL>IMULTBU)
557 if(!(IEL<200))std::cout <<
"5556:" << IEL << RHAZ << IMULTBU << std::endl;
558 if(IEL<0)std::cout <<
"5556:"<< IEL << RHAZ << IMULTBU << std::endl;
561 }
else if(IEL>IMULTBU){
570 }
else if(IEL>IMULTBU){
571 ATEMP =
dint(zprf + N_Breakup);
573 if((ATEMP - N_Breakup)<1.0 && N_Breakup<1.0){
585 aprf =
dint(zprf + N_Breakup);
587 NBU_DIFF = NBU_DIFF -
ISIGN(1,NBU_DIFF);
591 for(
G4int i=0;i<IMULTBU;i++)
595 if(ZBU_DIFF != 0 && NBU_DIFF == 0){
596 while(ZBU_DIFF > 0 || ZBU_DIFF < 0){
602 std::cout <<
"WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED." << std::endl;
606 if(IMEM_BU[IEL]==1)
goto mult5557;
608 if(IPROBA>IMULTBU+1 && ZBU_DIFF>0){
609 std::cout <<
"###',IPROBA,IMULTBU,NBU_DIFF,ZBU_DIFF,T_freeze_out" << std::endl;
617 N_Breakup = aprf - zprf;
619 aprf =
dint(zprf + N_Breakup);
624 if(!(IEL<200))std::cout <<
"5557:" << IEL << RHAZ << IMULTBU << std::endl;
625 if(IEL<0)std::cout <<
"5557:"<< IEL << RHAZ << IMULTBU << std::endl;
629 }
else if(IEL>IMULTBU){
630 N_Breakup =
dint(aprf - zprf);
633 ATEMP =
dint(ZTEMP + N_Breakup);
638 if((ATEMP-ZTEMP)<0.0){
642 if((ATEMP-ZTEMP)<1.0 && ZTEMP<1.0){
652 aprf =
dint(ZTEMP + N_Breakup);
655 ZBU_DIFF = ZBU_DIFF -
ISIGN(1,ZBU_DIFF);
665 for(
G4int i =0;i<IMULTBU;i++){
702 for(
G4int i = 0;i<IMULTBU;i++){
703 ABU_SUM = ABU_SUM +
BU_TAB[i][1];
704 ZBU_SUM = ZBU_SUM + BU_TAB[i][0];
706 ABU_DIFF =
idnint(ABU_SUM-AAINCL);
707 ZBU_DIFF =
idnint(ZBU_SUM-ZAINCL);
709 if(ABU_DIFF!=0 || ZBU_DIFF!=0)
710 std::cout <<
"Problem of mass in BU " << ABU_DIFF <<
" " << ZBU_DIFF << std::endl;
718 AMOMENT(AAINCL,aprf,1,&PXPRFP,&PYPRFP,&PZPRFP);
719 PPRFP = std::sqrt(PXPRFP*PXPRFP + PYPRFP*PYPRFP + PZPRFP*PZPRFP);
722 ETOT_PRF = std::sqrt(PPRFP*PPRFP + (aprf * amu)*(aprf * amu));
723 VX_PREF = C * PXPRFP / ETOT_PRF;
724 VY_PREF = C * PYPRFP / ETOT_PRF;
725 VZ_PREF = C * PZPRFP / ETOT_PRF;
728 tke_bu(zprf,aprf,ZAINCL,AAINCL,&VX1_BU,&VY1_BU,&VZ1_BU);
736 VX_PREF,VY_PREF,VZ_PREF,
737 &VXOUT,&VYOUT,&VZOUT);
744 VBU2 = VX_PREF*VX_PREF + VY_PREF*VY_PREF + VZ_PREF*VZ_PREF;
745 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C*C));
746 ETOT_PRF = aprf * amu / GAMMA_REL;
747 PXPRFP = ETOT_PRF * VX_PREF /
C;
748 PYPRFP = ETOT_PRF * VY_PREF /
C;
749 PZPRFP = ETOT_PRF * VZ_PREF /
C;
763 for(I_Breakup=0;I_Breakup<IMULTBU;I_Breakup++){
766 Eexc_BU_SUM = Eexc_BU_SUM +
BU_TAB[I_Breakup][2];
768 AMOMENT(AAINCL,BU_TAB[I_Breakup][1],1,&PX_BU,&PY_BU,&PZ_BU);
769 P_BU = std::sqrt(PX_BU*PX_BU + PY_BU*PY_BU + PZ_BU*PZ_BU);
772 ETOT_BU = std::sqrt(P_BU*P_BU + (BU_TAB[I_Breakup][1]*amu)*(BU_TAB[I_Breakup][1]*amu));
773 BU_TAB[I_Breakup][4] = C * PX_BU / ETOT_BU;
774 BU_TAB[I_Breakup][5] = C * PY_BU / ETOT_BU;
775 BU_TAB[I_Breakup][6] = C * PZ_BU / ETOT_BU;
777 tke_bu(BU_TAB[I_Breakup][0],BU_TAB[I_Breakup][1],ZAINCL,AAINCL,&VX2_BU,&VY2_BU,&VZ2_BU);
784 BU_TAB[I_Breakup][4],BU_TAB[I_Breakup][5],BU_TAB[I_Breakup][6],
785 &VXOUT,&VYOUT,&VZOUT);
787 BU_TAB[I_Breakup][4] = VXOUT;
788 BU_TAB[I_Breakup][5] = VYOUT;
789 BU_TAB[I_Breakup][6] = VZOUT;
792 VBU2 = BU_TAB[I_Breakup][4]*BU_TAB[I_Breakup][4] +
793 BU_TAB[I_Breakup][5]*BU_TAB[I_Breakup][5] +
794 BU_TAB[I_Breakup][6]*BU_TAB[I_Breakup][6];
795 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C*C));
796 ETOT_BU = BU_TAB[I_Breakup][1]*amu/GAMMA_REL;
797 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
798 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
799 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
801 PX_BU_SUM = PX_BU_SUM + PX_BU;
802 PY_BU_SUM = PY_BU_SUM + PY_BU;
803 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
808 P_BU_SUM = std::sqrt(PX_BU_SUM*PX_BU_SUM + PY_BU_SUM*PY_BU_SUM +
809 PZ_BU_SUM*PZ_BU_SUM);
812 ETOT_SUM = std::sqrt(P_BU_SUM*P_BU_SUM +
813 (AAINCL * amu)*(AAINCL * amu));
815 VX_BU_SUM = C * PX_BU_SUM / ETOT_SUM;
816 VY_BU_SUM = C * PY_BU_SUM / ETOT_SUM;
817 VZ_BU_SUM = C * PZ_BU_SUM / ETOT_SUM;
825 VX_PREF,VY_PREF,VZ_PREF,
826 &VXOUT,&VYOUT,&VZOUT);
832 VBU2 = VX_PREF*VX_PREF + VY_PREF*VY_PREF + VZ_PREF*VZ_PREF;
833 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C*C));
834 ETOT_PRF = aprf * amu / GAMMA_REL;
835 PXPRFP = ETOT_PRF * VX_PREF /
C;
836 PYPRFP = ETOT_PRF * VY_PREF /
C;
837 PZPRFP = ETOT_PRF * VZ_PREF /
C;
848 EKIN_BU = aprf * amu / GAMMA_REL - aprf *
amu;
850 for(I_Breakup=0;I_Breakup<IMULTBU;I_Breakup++){
858 &VXOUT,&VYOUT,&VZOUT);
860 BU_TAB[I_Breakup][4] = VXOUT;
861 BU_TAB[I_Breakup][5] = VYOUT;
862 BU_TAB[I_Breakup][6] = VZOUT;
867 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C*C));
869 ETOT_BU = BU_TAB[I_Breakup][1]*amu/GAMMA_REL;
871 EKIN_BU = EKIN_BU + BU_TAB[I_Breakup][1] * amu /
872 GAMMA_REL - BU_TAB[I_Breakup][1] *
amu;
874 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
875 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
876 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
877 E_tot_BU = E_tot_BU + ETOT_BU;
879 PX_BU_SUM = PX_BU_SUM + PX_BU;
880 PY_BU_SUM = PY_BU_SUM + PY_BU;
881 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
888 P_BU_SUM = std::sqrt(PX_BU_SUM*PX_BU_SUM + PY_BU_SUM*PY_BU_SUM +
889 PZ_BU_SUM*PZ_BU_SUM);
892 ETOT_SUM = std::sqrt(P_BU_SUM*P_BU_SUM +
893 (AAINCL * amu)*(AAINCL * amu));
895 VX_BU_SUM = C * PX_BU_SUM / ETOT_SUM;
896 VY_BU_SUM = C * PY_BU_SUM / ETOT_SUM;
897 VZ_BU_SUM = C * PZ_BU_SUM / ETOT_SUM;
905 VX_PREF,VY_PREF,VZ_PREF,
906 &VXOUT,&VYOUT,&VZOUT);
912 VBU2 = VX_PREF*VX_PREF + VY_PREF*VY_PREF + VZ_PREF*VZ_PREF;
913 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C*C));
914 ETOT_PRF = aprf * amu / GAMMA_REL;
915 PXPRFP = ETOT_PRF * VX_PREF /
C;
916 PYPRFP = ETOT_PRF * VY_PREF /
C;
917 PZPRFP = ETOT_PRF * VZ_PREF /
C;
928 EKIN_BU = aprf * amu / GAMMA_REL - aprf *
amu;
930 for(I_Breakup=0;I_Breakup<IMULTBU;I_Breakup++){
938 &VXOUT,&VYOUT,&VZOUT);
940 BU_TAB[I_Breakup][4] = VXOUT;
941 BU_TAB[I_Breakup][5] = VYOUT;
942 BU_TAB[I_Breakup][6] = VZOUT;
947 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C*C));
949 ETOT_BU = BU_TAB[I_Breakup][1]*amu/GAMMA_REL;
951 EKIN_BU = EKIN_BU + BU_TAB[I_Breakup][1] * amu /
952 GAMMA_REL - BU_TAB[I_Breakup][1] *
amu;
954 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
955 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
956 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
957 E_tot_BU = E_tot_BU + ETOT_BU;
959 PX_BU_SUM = PX_BU_SUM + PX_BU;
960 PY_BU_SUM = PY_BU_SUM + PY_BU;
961 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
969 for(
G4int i=0;i<IMULTBU;i++){
973 &VP1X,&VP1Y,&VP1Z,BU_TAB_TEMP,&ILOOP);
984 for(
int IJ=0;IJ<ILOOP;IJ++){
985 BU_TAB[IMULTBU+INEWLOOP+IJ][0] = BU_TAB_TEMP[IJ][0];
986 BU_TAB[IMULTBU+INEWLOOP+IJ][1] = BU_TAB_TEMP[IJ][1];
987 BU_TAB[IMULTBU+INEWLOOP+IJ][4] = BU_TAB_TEMP[IJ][2];
988 BU_TAB[IMULTBU+INEWLOOP+IJ][5] = BU_TAB_TEMP[IJ][3];
989 BU_TAB[IMULTBU+INEWLOOP+IJ][6] = BU_TAB_TEMP[IJ][4];
990 BU_TAB[IMULTBU+INEWLOOP+IJ][2] = 0.0;
991 BU_TAB[IMULTBU+INEWLOOP+IJ][3] = 0.0;
994 INEWLOOP = INEWLOOP + ILOOP;
1001 IMULTBU = IMULTBU + INEWLOOP;
1014 for(
G4int i=0;i<IMULTBU;i++){
1018 Nblamb =
new G4int[IMULTBU];
1019 for(
G4int i=0;i<IMULTBU;i++)Nblamb[i] = 0;
1020 for(
G4int j=0;j<NbLam0;){
1021 G4double probtotal = (aprf - zprf)/sumN;
1024 if(ran <= probtotal){
1028 for(
G4int i=0;i<IMULTBU;i++){
1030 if(probtotal < ran && ran <= probtotal+problamb[i]){
1031 Nblamb[i] = Nblamb[i] + 1;
1034 probtotal = probtotal + problamb[i];
1040 for(
G4int i=0;i<IMULTBU;i++){
1045 G4int nbl = Nblamb[i];
1046 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,&nbl);
1052 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1053 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1054 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1055 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1062 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1063 &VXOUT,&VYOUT,&VZOUT);
1064 EV_TAB[IJ+IEV_TAB][2] = VXOUT;
1065 EV_TAB[IJ+IEV_TAB][3] = VYOUT;
1066 EV_TAB[IJ+IEV_TAB][4] = VZOUT;
1068 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1078 &VXOUT,&VYOUT,&VZOUT);
1098 G4double EkinR1 = TKEIMFBU * AIMFBU / (AFBU+AIMFBU);
1099 G4double EkinR2 = TKEIMFBU * AFBU / (AFBU+AIMFBU);
1100 G4double V1 = std::sqrt(EkinR1/AFBU) * 1.3887;
1101 G4double V2 = std::sqrt(EkinR2/AIMFBU) * 1.3887;
1103 G4double VPERP1 = std::sqrt(V1*V1 - VZ1_IMF*VZ1_IMF);
1105 G4double VX1_IMF = VPERP1 * std::sin(ALPHA1);
1106 G4double VY1_IMF = VPERP1 * std::cos(ALPHA1);
1107 G4double VX2_IMF = - VX1_IMF / V1 * V2;
1108 G4double VY2_IMF = - VY1_IMF / V1 * V2;
1109 G4double VZ2_IMF = - VZ1_IMF / V1 * V2;
1111 G4double EEIMFP = EEBU * AFBU /(AFBU + AIMFBU);
1112 G4double EEIMF = EEBU * AIMFBU /(AFBU + AIMFBU);
1115 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.));
1117 G4double JPRFHEAVY =
BU_TAB[i][9] * 0.4 * 931.49 * 1.16*1.16 * std::pow(AFBU,5.0/3.0) / IINERTTOT;
1118 G4double JPRFLIGHT =
BU_TAB[i][9] * 0.4 * 931.49 * 1.16*1.16 * std::pow(AIMFBU,5.0/3.0) / IINERTTOT;
1127 &VXOUT,&VYOUT,&VZOUT);
1132 G4double vx1ev_imf=0., vy1ev_imf=0., vz1ev_imf=0., zdummy=0., adummy=0., tkedummy=0.,jprf1=0.;
1137 G4double pbH = (AFBU-ZFBU) / (AFBU-ZFBU+AIMFBU-ZIMFBU);
1138 for(
G4int j=0;j<nbl;j++){
1146 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,&NbLamH);
1148 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1149 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1150 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1151 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1158 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1159 &VXOUT,&VYOUT,&VZOUT);
1160 EV_TAB[IJ+IEV_TAB][2] = VXOUT;
1161 EV_TAB[IJ+IEV_TAB][3] = VYOUT;
1162 EV_TAB[IJ+IEV_TAB][4] = VZOUT;
1164 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1175 &VXOUT,&VYOUT,&VZOUT);
1185 G4double zffimf, affimf,zdummy1, adummy1, tkedummy1, jprf2, vx2ev_imf, vy2ev_imf, vz2ev_imf;
1187 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,&NbLamimf);
1189 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1190 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1191 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1192 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1199 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1200 &VXOUT,&VYOUT,&VZOUT);
1203 &VX2OUT,&VY2OUT,&VZ2OUT);
1204 EV_TAB[IJ+IEV_TAB][2] = VX2OUT;
1205 EV_TAB[IJ+IEV_TAB][3] = VY2OUT;
1206 EV_TAB[IJ+IEV_TAB][4] = VZ2OUT;
1208 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1216 BU_TAB[IMULTBU+ILOOPBU][11]= NbLamimf;
1220 &VXOUT,&VYOUT,&VZOUT);
1223 &VX2OUT,&VY2OUT,&VZ2OUT);
1224 BU_TAB[IMULTBU+ILOOPBU][4] = VX2OUT;
1225 BU_TAB[IMULTBU+ILOOPBU][5] = VY2OUT;
1226 BU_TAB[IMULTBU+ILOOPBU][6] = VZ2OUT;
1227 ILOOPBU = ILOOPBU + 1;
1240 BU_TAB[i][11]= Nblamb[i];
1244 IMULTBU = IMULTBU + ILOOPBU;
1252 for(
G4int i=0;i<IMULTBU;i++){
1253 ABU_SUM = ABU_SUM +
BU_TAB[i][8];
1254 ZBU_SUM = ZBU_SUM + BU_TAB[i][7];
1256 BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6],
1257 &VP1X,&VP1Y,&VP1Z,BU_TAB_TEMP1,&ILOOP);
1265 ABU_SUM = ABU_SUM +
G4double(afpnew) - BU_TAB[i][8];
1266 ZBU_SUM = ZBU_SUM +
G4double(zfpnew) - BU_TAB[i][7];
1269 BU_TAB[i][4] = VP1X;
1270 BU_TAB[i][5] = VP1Y;
1271 BU_TAB[i][6] = VP1Z;
1274 for(
G4int IJ=0;IJ<ILOOP;IJ++){
1275 BU_TAB[IMULTBU+INEWLOOP+IJ][7] = BU_TAB_TEMP1[IJ][0];
1276 BU_TAB[IMULTBU+INEWLOOP+IJ][8] = BU_TAB_TEMP1[IJ][1];
1277 BU_TAB[IMULTBU+INEWLOOP+IJ][4] = BU_TAB_TEMP1[IJ][2];
1278 BU_TAB[IMULTBU+INEWLOOP+IJ][5] = BU_TAB_TEMP1[IJ][3];
1279 BU_TAB[IMULTBU+INEWLOOP+IJ][6] = BU_TAB_TEMP1[IJ][4];
1280 BU_TAB[IMULTBU+INEWLOOP+IJ][2] = 0.0;
1281 BU_TAB[IMULTBU+INEWLOOP+IJ][3] = 0.0;
1282 BU_TAB[IMULTBU+INEWLOOP+IJ][0] = BU_TAB[i][0];
1283 BU_TAB[IMULTBU+INEWLOOP+IJ][1] = BU_TAB[i][1];
1284 BU_TAB[IMULTBU+INEWLOOP+IJ][11] = BU_TAB[i][11];
1285 ABU_SUM = ABU_SUM + BU_TAB[IMULTBU+INEWLOOP+IJ][8];
1286 ZBU_SUM = ZBU_SUM + BU_TAB[IMULTBU+INEWLOOP+IJ][7];
1289 INEWLOOP = INEWLOOP + ILOOP;
1294 IMULTBU = IMULTBU + INEWLOOP;
1298 VX_PREF,VY_PREF,VZ_PREF,
1299 &VXOUT,&VYOUT,&VZOUT);
1304 for(
G4int i=0;i<IMULTBU;i++){
1307 &VXOUT,&VYOUT,&VZOUT);
1312 for(
G4int i=0;i<IEV_TAB;i++){
1315 &VXOUT,&VYOUT,&VZOUT);
1322 if(IMULTBU>200)std::cout <<
"IMULTBU>200 " << IMULTBU << std::endl;
1353 if(zprfp<=2 && zprfp<aprfp){
1371 if(zprfp<=2 && zprfp==aprfp){
1373 VX_PREF, VY_PREF, VZ_PREF,
1374 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
1381 for(
G4int IJ = 0; IJ<6; IJ++)
1382 EV_TAB[
I+IEV_TAB][IJ] = EV_TAB_TEMP[
I][IJ];
1384 IEV_TAB = IEV_TAB + ILOOP;
1402 VX_PREF, VY_PREF, VZ_PREF,
1403 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
1410 for(
G4int IJ = 0; IJ<6; IJ++)
1411 EV_TAB[
I+IEV_TAB][IJ] = EV_TAB_TEMP[
I][IJ];
1413 IEV_TAB = IEV_TAB + ILOOP;
1428 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,&NbLam0);
1430 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1431 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1432 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1433 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1440 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1441 &VXOUT,&VYOUT,&VZOUT);
1442 EV_TAB[IJ+IEV_TAB][2] = VXOUT;
1443 EV_TAB[IJ+IEV_TAB][3] = VYOUT;
1444 EV_TAB[IJ+IEV_TAB][4] = VZOUT;
1446 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1452 vx_eva,vy_eva,vz_eva,
1453 &VXOUT,&VYOUT,&VZOUT);
1458 if(ff == 0 && fimf == 0){
1470 VFP1_CM[0] = V_CM[0];
1471 VFP1_CM[1] = V_CM[1];
1472 VFP1_CM[2] = V_CM[2];
1473 for(
G4int j=0;j<3;j++){
1479 if(ff == 1 && fimf == 0) ftype = 1;
1480 if(ff == 0 && fimf == 1) ftype = 2;
1493 G4int IEV_TAB_FIS = 0,imode=0;
1495 G4double vx1_fission=0.,vy1_fission=0.,vz1_fission=0.;
1496 G4double vx2_fission=0.,vy2_fission=0.,vz2_fission=0.;
1497 G4double vx_eva_sc=0.,vy_eva_sc=0.,vz_eva_sc=0.;
1500 &vx1_fission,&vy1_fission,&vz1_fission,
1501 &vx2_fission,&vy2_fission,&vz2_fission,
1502 &ZFP1,&AFP1,&SFP1,&ZFP2,&AFP2,&SFP2,&imode,
1503 &vx_eva_sc,&vy_eva_sc,&vz_eva_sc,EV_TEMP,&IEV_TAB_FIS,&NbLam0);
1505 for(
G4int IJ = 0; IJ< IEV_TAB_FIS;IJ++){
1506 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1507 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1508 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1515 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1516 &VXOUT,&VYOUT,&VZOUT);
1517 EV_TAB[IJ+IEV_TAB][2] = VXOUT;
1518 EV_TAB[IJ+IEV_TAB][3] = VYOUT;
1519 EV_TAB[IJ+IEV_TAB][4] = VZOUT;
1521 IEV_TAB = IEV_TAB + IEV_TAB_FIS;
1537 V_CM[0],V_CM[1],V_CM[2],
1538 &VXOUT,&VYOUT,&VZOUT);
1541 &VX2OUT,&VY2OUT,&VZ2OUT);
1542 VFP1_CM[0] = VX2OUT;
1543 VFP1_CM[1] = VY2OUT;
1544 VFP1_CM[2] = VZ2OUT;
1551 V_CM[0],V_CM[1],V_CM[2],
1552 &VXOUT,&VYOUT,&VZOUT);
1555 &VX2OUT,&VY2OUT,&VZ2OUT);
1556 VFP2_CM[0] = VX2OUT;
1557 VFP2_CM[1] = VY2OUT;
1558 VFP2_CM[2] = VZ2OUT;
1562 }
else if(ftype == 2){
1571 G4double pbH = (af-zf) / (af-zf+aimf-zimf);
1573 for(
G4int i=0;i<NbLam0;i++){
1582 G4double EkinR1 = tkeimf * aimf / (af+aimf);
1583 G4double EkinR2 = tkeimf * af / (af+aimf);
1585 G4double V2 = std::sqrt(EkinR2/aimf) * 1.3887;
1587 G4double VPERP1 = std::sqrt(V1*V1 - VZ1_IMF*VZ1_IMF);
1589 G4double VX1_IMF = VPERP1 * std::sin(ALPHA1);
1590 G4double VY1_IMF = VPERP1 * std::cos(ALPHA1);
1591 G4double VX2_IMF = - VX1_IMF / V1 * V2;
1592 G4double VY2_IMF = - VY1_IMF / V1 * V2;
1593 G4double VZ2_IMF = - VZ1_IMF / V1 * V2;
1595 G4double EEIMFP = ee * af /(af + aimf);
1596 G4double EEIMF = ee * aimf /(af + aimf);
1599 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.));
1601 G4double JPRFHEAVY = jprf0 * 0.4 * 931.49 * 1.16*1.16 * std::pow(af,5.0/3.0) / IINERTTOT;
1602 G4double JPRFLIGHT = jprf0 * 0.4 * 931.49 * 1.16*1.16 * std::pow(aimf,5.0/3.0) / IINERTTOT;
1603 if(af<2.0) std::cout <<
"RN117-4,AF,ZF,EE,JPRFheavy" << std::endl;
1605 G4double vx1ev_imf=0., vy1ev_imf=0., vz1ev_imf=0., zdummy=0., adummy=0., tkedummy=0.,jprf1=0.;
1607 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,&NbLamH);
1609 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1610 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1611 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1612 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1619 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1620 &VXOUT,&VYOUT,&VZOUT);
1623 &VX2OUT,&VY2OUT,&VZ2OUT);
1624 EV_TAB[IJ+IEV_TAB][2] = VX2OUT;
1625 EV_TAB[IJ+IEV_TAB][3] = VY2OUT;
1626 EV_TAB[IJ+IEV_TAB][4] = VZ2OUT;
1628 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1637 G4double zffimf, affimf,zdummy1=0., adummy1=0., tkedummy1=0.,jprf2,vx2ev_imf,vy2ev_imf,
1640 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,&NbLamimf);
1642 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1643 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1644 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1645 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1652 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1653 &VXOUT,&VYOUT,&VZOUT);
1656 &VX2OUT,&VY2OUT,&VZ2OUT);
1657 EV_TAB[IJ+IEV_TAB][2] = VX2OUT;
1658 EV_TAB[IJ+IEV_TAB][3] = VY2OUT;
1659 EV_TAB[IJ+IEV_TAB][4] = VZ2OUT;
1661 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1675 V_CM[0],V_CM[1],V_CM[2],
1676 &VXOUT,&VYOUT,&VZOUT);
1679 &VX2OUT,&VY2OUT,&VZ2OUT);
1680 VIMF_CM[0] = VX2OUT;
1681 VIMF_CM[1] = VY2OUT;
1682 VIMF_CM[2] = VZ2OUT;
1688 V_CM[0],V_CM[1],V_CM[2],
1689 &VXOUT,&VYOUT,&VZOUT);
1692 &VX2OUT,&VY2OUT,&VZ2OUT);
1693 VFP1_CM[0] = VX2OUT;
1694 VFP1_CM[1] = VY2OUT;
1695 VFP1_CM[2] = VZ2OUT;
1697 if(FF11==0 && FIMF11==0){
1713 }
else if(FF11==1 && FIMF11==0){
1727 G4int IEV_TAB_FIS = 0,imode=0;
1729 G4double vx1_fission=0.,vy1_fission=0.,vz1_fission=0.;
1730 G4double vx2_fission=0.,vy2_fission=0.,vz2_fission=0.;
1731 G4double vx_eva_sc=0.,vy_eva_sc=0.,vz_eva_sc=0.;
1734 &vx1_fission,&vy1_fission,&vz1_fission,
1735 &vx2_fission,&vy2_fission,&vz2_fission,
1736 &ZFP1,&AFP1,&SFP1,&ZFP2,&AFP2,&SFP2,&imode,
1737 &vx_eva_sc,&vy_eva_sc,&vz_eva_sc,EV_TEMP,&IEV_TAB_FIS,&NbLamH);
1739 for(
int IJ = 0; IJ< IEV_TAB_FIS;IJ++){
1740 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1741 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1742 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1749 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1750 &VXOUT,&VYOUT,&VZOUT);
1751 EV_TAB[IJ+IEV_TAB][2] = VXOUT;
1752 EV_TAB[IJ+IEV_TAB][3] = VYOUT;
1753 EV_TAB[IJ+IEV_TAB][4] = VZOUT;
1755 IEV_TAB = IEV_TAB + IEV_TAB_FIS;
1768 V_CM[0],V_CM[1],V_CM[2],
1769 &VXOUT,&VYOUT,&VZOUT);
1772 &VX2OUT,&VY2OUT,&VZ2OUT);
1774 VX2OUT,VY2OUT,VZ2OUT,
1775 &VXOUT,&VYOUT,&VZOUT);
1778 &VX2OUT,&VY2OUT,&VZ2OUT);
1779 VFP1_CM[0] = VX2OUT;
1780 VFP1_CM[1] = VY2OUT;
1781 VFP1_CM[2] = VZ2OUT;
1791 V_CM[0],V_CM[1],V_CM[2],
1792 &VXOUT,&VYOUT,&VZOUT);
1795 &VX2OUT,&VY2OUT,&VZ2OUT);
1797 VX2OUT,VY2OUT,VZ2OUT,
1798 &VXOUT,&VYOUT,&VZOUT);
1801 &VX2OUT,&VY2OUT,&VZ2OUT);
1802 VFP2_CM[0] = VX2OUT;
1803 VFP2_CM[1] = VY2OUT;
1804 VFP2_CM[2] = VZ2OUT;
1806 }
else if(FF11==0 && FIMF11==1){
1823 G4double pbH1 = (af-zf) / (af-zf+aimf-zimf);
1824 for(
G4int i=0;i<NbLamH;i++){
1833 EkinR1 = tkeimf * aimf / (af+aimf);
1834 EkinR2 = tkeimf * af / (af+aimf);
1835 V1 = std::sqrt(EkinR1/af) * 1.3887;
1836 V2 = std::sqrt(EkinR2/aimf) * 1.3887;
1838 VPERP1 = std::sqrt(V1*V1 - VZ1_IMFS*VZ1_IMFS);
1840 G4double VX1_IMFS = VPERP1 * std::sin(ALPHA1);
1841 G4double VY1_IMFS = VPERP1 * std::cos(ALPHA1);
1842 G4double VX2_IMFS = - VX1_IMFS / V1 * V2;
1843 G4double VY2_IMFS = - VY1_IMFS / V1 * V2;
1844 G4double VZ2_IMFS = - VZ1_IMFS / V1 * V2;
1846 EEIMFP = ee * af /(af + aimf);
1847 EEIMF = ee * aimf /(af + aimf);
1850 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.));
1852 JPRFHEAVY = jprf1 * 0.4 * 931.49 * 1.16*1.16 * std::pow(af,5.0/3.0) / IINERTTOT;
1853 JPRFLIGHT = jprf1 * 0.4 * 931.49 * 1.16*1.16 * std::pow(aimf,5.0/3.0) / IINERTTOT;
1855 G4double zffs=0.,affs=0.,vx1ev_imfs=0.,vy1ev_imfs=0.,vz1ev_imfs=0.,jprf3=0.;
1857 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,&NbLamH1);
1859 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1860 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1861 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1862 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1869 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1870 &VXOUT,&VYOUT,&VZOUT);
1873 &VX2OUT,&VY2OUT,&VZ2OUT);
1874 EV_TAB[IJ+IEV_TAB][2] = VX2OUT;
1875 EV_TAB[IJ+IEV_TAB][3] = VY2OUT;
1876 EV_TAB[IJ+IEV_TAB][4] = VZ2OUT;
1878 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1887 G4double zffimfs=0.,affimfs=0.,vx2ev_imfs=0.,vy2ev_imfs=0.,vz2ev_imfs=0.,jprf4=0.;
1889 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,&NbLamimf1);
1891 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1892 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1893 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1894 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1901 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1902 &VXOUT,&VYOUT,&VZOUT);
1905 &VX2OUT,&VY2OUT,&VZ2OUT);
1906 EV_TAB[IJ+IEV_TAB][2] = VX2OUT;
1907 EV_TAB[IJ+IEV_TAB][3] = VY2OUT;
1908 EV_TAB[IJ+IEV_TAB][4] = VZ2OUT;
1910 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1925 V_CM[0],V_CM[1],V_CM[2],
1926 &VXOUT,&VYOUT,&VZOUT);
1929 &VX2OUT,&VY2OUT,&VZ2OUT);
1931 VX2OUT,VY2OUT,VZ2OUT,
1932 &VXOUT,&VYOUT,&VZOUT);
1935 &VX2OUT,&VY2OUT,&VZ2OUT);
1936 VFP1_CM[0] = VX2OUT;
1937 VFP1_CM[1] = VY2OUT;
1938 VFP1_CM[2] = VZ2OUT;
1946 V_CM[0],V_CM[1],V_CM[2],
1947 &VXOUT,&VYOUT,&VZOUT);
1950 &VX2OUT,&VY2OUT,&VZ2OUT);
1952 VX2OUT,VY2OUT,VZ2OUT,
1953 &VXOUT,&VYOUT,&VZOUT);
1956 &VX2OUT,&VY2OUT,&VZ2OUT);
1957 VFP2_CM[0] = VX2OUT;
1958 VFP2_CM[1] = VY2OUT;
1959 VFP2_CM[2] = VZ2OUT;
1964 if(ftype!=1 && ftype!=21){
1970 VFP1_CM[0],VFP1_CM[1],VFP1_CM[2],
1971 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
1980 for(
G4int IJ = 0; IJ<5; IJ++)
1981 EV_TAB[
I+IEV_TAB][IJ] = EV_TAB_TEMP[
I][IJ];
1983 IEV_TAB = IEV_TAB + ILOOP;
1990 VIMF_CM[0],VIMF_CM[1],VIMF_CM[2],
1991 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
2000 for(
G4int IJ = 0; IJ<5; IJ++)
2001 EV_TAB[
I+IEV_TAB][IJ] = EV_TAB_TEMP[
I][IJ];
2003 IEV_TAB = IEV_TAB + ILOOP;
2010 VFP2_CM[0],VFP2_CM[1],VFP2_CM[2],
2011 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
2020 for(
G4int IJ = 0; IJ<5; IJ++)
2021 EV_TAB[
I+IEV_TAB][IJ] = EV_TAB_TEMP[
I][IJ];
2023 IEV_TAB = IEV_TAB + ILOOP;
2031 if(ftype==1 || ftype==21){
2036 VFP1_CM[0],VFP1_CM[1],VFP1_CM[2],
2037 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
2046 for(
G4int IJ = 0; IJ<5; IJ++)
2047 EV_TAB[
I+IEV_TAB][IJ] = EV_TAB_TEMP[
I][IJ];
2049 IEV_TAB = IEV_TAB + ILOOP;
2055 VFP2_CM[0],VFP2_CM[1],VFP2_CM[2],
2056 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
2065 for(
G4int IJ = 0; IJ<5; IJ++)
2066 EV_TAB[
I+IEV_TAB][IJ] = EV_TAB_TEMP[
I][IJ];
2068 IEV_TAB = IEV_TAB + ILOOP;
2075 VIMF_CM[0],VIMF_CM[1],VIMF_CM[2],
2076 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
2085 for(
G4int IJ = 0; IJ<5; IJ++)
2086 EV_TAB[
I+IEV_TAB][IJ] = EV_TAB_TEMP[
I][IJ];
2088 IEV_TAB = IEV_TAB + ILOOP;
2094 if((ftype == 1 || ftype == 21) && (AFP2<=0 || AFP1<=0 || ZFP2<=0 || ZFP1<=0)){
2095 std::cout <<
"ZFP1:" << ZFP1 << std::endl;
2096 std::cout <<
"AFP1:" << AFP1 << std::endl;
2097 std::cout <<
"ZFP2:" << ZFP2 << std::endl;
2098 std::cout <<
"AFP2:" << AFP2 << std::endl;
2102 EV_TAB[IEV_TAB][0] = ZFP1;
2103 EV_TAB[IEV_TAB][1] = AFP1;
2104 EV_TAB[IEV_TAB][5] = SFP1;
2105 EV_TAB[IEV_TAB][2] = VFP1_CM[0];
2106 EV_TAB[IEV_TAB][3] = VFP1_CM[1];
2107 EV_TAB[IEV_TAB][4] = VFP1_CM[2];
2108 IEV_TAB = IEV_TAB + 1;
2111 EV_TAB[IEV_TAB][0] = ZFP2;
2112 EV_TAB[IEV_TAB][1] = AFP2;
2113 EV_TAB[IEV_TAB][5] = SFP2;
2114 EV_TAB[IEV_TAB][2] = VFP2_CM[0];
2115 EV_TAB[IEV_TAB][3] = VFP2_CM[1];
2116 EV_TAB[IEV_TAB][4] = VFP2_CM[2];
2117 IEV_TAB = IEV_TAB + 1;
2121 EV_TAB[IEV_TAB][0] = ZFPIMF;
2122 EV_TAB[IEV_TAB][1] = AFPIMF;
2123 EV_TAB[IEV_TAB][5] = SFPIMF;
2124 EV_TAB[IEV_TAB][2] = VIMF_CM[0];
2125 EV_TAB[IEV_TAB][3] = VIMF_CM[1];
2126 EV_TAB[IEV_TAB][4] = VIMF_CM[2];
2127 IEV_TAB = IEV_TAB + 1;
2191 #ifdef ABLAXX_IN_GEANT4_MODE
2196 if(dataInterface->
readData() ==
true) {
2234 for(
G4int i=1;i<13;i++){
2235 for(
G4int j=1;j<154;j++){
2245 mfrldm[j][i] = MP*i+MN*j+
eflmac(i+j,i,1,0);
2250 for(
G4int i=1;i<13;i++){
2251 for(
G4int j=1;j<154;j++){
2271 e0 = 0.285+11.17*std::pow(j+i,-0.464) -0.390-0.00058*(j+i);
2277 e0 = 22.34*std::pow(j+i,-0.464)-0.235;
2284 if((j==i)&&
mod(j,2)==1&&
mod(i,2)==1){
2285 e0 = e0 - 30.0*(1.0/
G4double(j+i));
2299 delete dataInterface;
2394 G4double xv = 0.0, xs = 0.0, xc = 0.0, xa = 0.0;
2396 if ((a <= 0.01) || (z < 0.01)) {
2401 xs = 17.23*std::pow(a,(2.0/3.0));
2404 xc = 0.7*z*(z-1.0)*std::pow((a-1.0),(-1.e0/3.e0));
2411 xa = 23.6*(std::pow((a-2.0*z),2)/
a);
2412 (*el) = xv+xs+xc+xa;
2440 if ( (a1 <= 0) || (z1 <= 0) || ((a1-z1) <= 0) ) {
2449 (*el) =
eflmac(a1,z1,0,refopt4);
2459 (*el) = (*el) + (12.552-0.1436*
z1);
2461 if(n1>145&&n1<=152){
2462 (*el) = (*el) + ((152.4-1.77*
z1)+(-0.972+0.0113*z1)*
n1);
2484 const G4int alpha2Size = 37;
2486 G4double alpha2[alpha2Size] = {0.0, 2.5464e0, 2.4944e0, 2.4410e0, 2.3915e0, 2.3482e0,
2487 2.3014e0, 2.2479e0, 2.1982e0, 2.1432e0, 2.0807e0, 2.0142e0, 1.9419e0,
2488 1.8714e0, 1.8010e0, 1.7272e0, 1.6473e0, 1.5601e0, 1.4526e0, 1.3164e0,
2489 1.1391e0, 0.9662e0, 0.8295e0, 0.7231e0, 0.6360e0, 0.5615e0, 0.4953e0,
2490 0.4354e0, 0.3799e0, 0.3274e0, 0.2779e0, 0.2298e0, 0.1827e0, 0.1373e0,
2491 0.0901e0, 0.0430e0, 0.0000e0};
2496 v = (x - 0.3)/
dx + 1.0;
2507 return(alpha2[index] + (alpha2[index+1] - alpha2[index]) /
dx * ( x - (0.3e0 +
dx*(index-1))));
2522 G4double aa = 0.0,
zz = 0.0, i = 0.0,z2a,C_S,
R,W,
G,G1,G2,A_CC;
2528 z2a=
zz*
zz/aa-ny*(1115.-939.+sn-slam)/(0.7053*std::pow(a,2./3.));
2532 fissilityResult = std::pow(
zz,2) / aa /50.8830e0 / (1.0e0 - 1.7826e0 * std::pow(i,2));
2537 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));
2542 fissilityResult = std::pow(
zz,2) / aa /(48.e0*(1.e0 - 17.22e0*std::pow(i,4)));
2547 C_S = 21.13 * (1.0 - 2.3*i*i);
2548 R = 1.16 * std::pow(aa,1.0/3.0);
2550 G1 = 1.0 - 15.0/8.0*W+21.0/8.0*W*W*W;
2551 G2 = 1.0 + 9.0/2.0*W + 7.0*W*W + 7.0/2.0*W*W*W;
2552 G = 1.0 - 5.0*W*W*(G1 - 3.0/4.0*G2*std::exp(-2.0/W));
2553 A_CC = 3.0/5.0 * 1.44 * G / 1.16;
2554 fissilityResult = z2a * A_CC/(2.0*C_S);
2557 if (fissilityResult > 1.0) {
2558 fissilityResult = 1.0;
2561 if (fissilityResult < 0.0) {
2562 fissilityResult = 0.0;
2565 return fissilityResult;
2568 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,
2569 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][6],
G4int *iev_tab_temp_par,
G4int *NbLam0_par)
2580 G4int fimf = (*fimf_par);
2582 G4int inttype = (*inttype_par);
2583 G4int inum = (*inum_par);
2584 G4int NbLam0 = (*NbLam0_par);
2676 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, problamb0 = 0.0, ptotl = 0.0,
e = 0.0, tcn = 0.0;
2677 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, slamb0 = 0.0;
2678 G4double ecn = 0.0, ecp = 0.0, ecd = 0.0, ect = 0.0,eche = 0.0,eca = 0.0, ecg = 0.0, eclamb0 = 0.0,
bp = 0.0, bd = 0.0, bt = 0.0, bhe = 0.0, ba = 0.0;
2679 G4double zimf= 0.0,aimf= 0.0,bimf= 0.0,sbimf= 0.0,timf= 0.0;
2682 G4double ctet1 = 0.0, stet1 = 0.0, phi1 = 0.0;
2686 G4int fgamma = 0, gammadecay = 0, flamb0decay=0;
2688 G4double jprfn=0.0, jprfp=0.0, jprfd=0.0, jprft=0.0, jprfhe=0.0, jprfa=0.0, jprflamb0 = 0.0;
2694 const G4double mu2 = 931.494*931.494;
2699 G4int IEV_TAB_TEMP=0;
2701 for(
G4int I1=0;I1<200;I1++)
2702 for(
G4int I2=0;I2<6;I2++)
2703 EV_TEMP[I1][I2] = 0.0;
2713 if(ee<0.|| zf<3.)
goto evapora100;
2714 direct(zf,af,ee,jprf,&probp,&probd,&probt,&probn,&probhe,&proba,&probg,&probimf,&probf,&problamb0,&ptotl,
2715 &sn,&sbp,&sbd,&sbt,&sbhe,&sba,&slamb0,
2716 &ecn,&ecp,&ecd,&ect,&eche,&eca,&ecg,&eclamb0,
2717 &
bp,&bd,&bt,&bhe,&ba,&
sp,&sd,&st,&she,&sa,&ef,&ts1,inttype,inum,itest,&sortie,&tcn,
2718 &jprfn, &jprfp, &jprfd, &jprft, &jprfhe, &jprfa, &jprflamb0, &tsum, NbLam0);
2722 if(ptotl==0.0)
goto evapora100;
2726 if(
e>1e30)std::cout <<
"ERROR AT THE EXIT OF EVAPORA,E>1.D30,AF="<< af <<
" ZF=" << zf << std::endl;
2733 pc = std::sqrt(std::pow((1.0 + (ecn)/9.3956
e2),2.) - 1.0) * 9.3956e2;
2740 else if(probp!=0.0){
2744 pc = std::sqrt(std::pow((1.0 + ecp/9.3827
e2),2.) - 1.0) * 9.3827e2;
2751 else if(probd!=0.0){
2755 pc = std::sqrt(std::pow((1.0 + ecd/1.875358
e3),2) - 1.0) * 1.875358e3;
2762 else if(probt!=0.0){
2766 pc = std::sqrt(std::pow((1.0 + ect/2.80828
e3),2) - 1.0) * 2.80828e3;
2773 else if(probhe!=0.0){
2776 epsiln = she + eche;
2777 pc = std::sqrt(std::pow((1.0 + eche/2.80826
e3),2) - 1.0) * 2.80826e3;
2784 else{
if(proba!=0.0){
2788 pc = std::sqrt(std::pow((1.0 + eca/3.72834
e3),2) - 1.0) * 3.72834e3;
2810 pc = std::sqrt(std::pow((1.0 + eca/3.72834
e3),2) - 1.0) * 3.72834e3;
2819 else if (
x < proba+probhe) {
2823 epsiln = she + eche;
2824 pc = std::sqrt(std::pow((1.0 + eche/2.80826
e3),2) - 1.0) * 2.80826e3;
2833 else if (
x < proba+probhe+probt) {
2838 pc = std::sqrt(std::pow((1.0 + ect/2.80828
e3),2) - 1.0) * 2.80828e3;
2847 else if (
x < proba+probhe+probt+probd) {
2852 pc = std::sqrt(std::pow((1.0 + ecd/1.875358
e3),2) - 1.0) * 1.875358e3;
2861 else if (
x < proba+probhe+probt+probd+probp) {
2866 pc = std::sqrt(std::pow((1.0 + ecp/9.3827
e2),2) - 1.0) * 9.3827e2;
2875 else if (
x < proba+probhe+probt+probd+probp+probn) {
2880 pc = std::sqrt(std::pow((1.0 + (ecn)/9.3956
e2),2.) - 1.0) * 9.3956e2;
2889 else if (
x < proba+probhe+probt+probd+probp+probn+problamb0) {
2893 epsiln = slamb0 + eclamb0;
2894 pc = std::sqrt(std::pow((1.0 + (eclamb0)/11.1568
e2),2.) - 1.0) * 11.1568e2;
2905 else if (
x < proba+probhe+probt+probd+probp+probn+problamb0+probg) {
2916 if(probp==0.0 && probn==0.0 && probd==0.0 && probt==0.0 && proba==0.0 && probhe==0.0 && problamb0==0.0 && probimf==0.0 && probf==0.0)fgamma = 1;
2920 else if (
x < proba+probhe+probt+probd+probp+probn+problamb0+probg+probimf) {
2927 imf(af,zf,tcn,ee,&zimf,&aimf,&bimf,&sbimf,&timf,jprf);
2929 if(iloop>100)std::cout <<
"Problem in EVAPORA: IMF called > 100 times" << std::endl;
2930 if(zimf>=(zf-2.0))
goto dir1973;
2936 if(zimf==0.0 || aimf==0.0 || sbimf>ee)std::cout <<
"warning: Look in EVAPORA CALL IMF" << std::endl;
2947 tkeimf=
min(2.0*timf,ee-sbimf);
2950 if(tkeimf<=0.0)
goto dir1235;
2951 if(tkeimf>(ee-sbimf) && timf>0.5)
goto dir1235;
2953 tkeimf = tkeimf + bimf;
2957 epsiln = (sbimf-bimf) + tkeimf;
2986 if (ee <= 0.01)ee = 0.01;
2988 if(gammadecay==1 && ee<(epsiln+0.010)){
2989 epsiln = ee - 0.010;
2994 std::cout <<
"***WARNING epsilon<0***" << std::endl;
3002 if (ee <= 0.01)ee = 0.01;
3003 mtota = mtota + malpha;
3010 if(amoins==2 && zmoins==0){twon=1;amoins=1;}
else{ twon=0;}
3014 if(ff==0 && fimf==0){
3017 EV_TEMP[IEV_TAB_TEMP][0] = 0.;
3018 EV_TEMP[IEV_TAB_TEMP][1] = -2;
3019 EV_TEMP[IEV_TAB_TEMP][5] = 1.;
3021 EV_TEMP[IEV_TAB_TEMP][0] = zmoins;
3022 EV_TEMP[IEV_TAB_TEMP][1] = amoins;
3023 EV_TEMP[IEV_TAB_TEMP][5] = 0.;
3026 ctet1 = 2.0*rnd - 1.0;
3027 stet1 = std::sqrt(1.0 - std::pow(ctet1,2));
3029 phi1 = rnd*2.0*3.141592654;
3030 G4double xcv = stet1*std::cos(phi1);
3031 G4double ycv = stet1*std::sin(phi1);
3036 G4double ETOT_LP = std::sqrt(pc*pc + amoins*amoins * mu2);
3037 if(flamb0decay==1)ETOT_LP = std::sqrt(pc*pc + 1115.683*1115.683);
3038 EV_TEMP[IEV_TAB_TEMP][2] = c * pc * xcv / ETOT_LP;
3039 EV_TEMP[IEV_TAB_TEMP][3] = c * pc * ycv / ETOT_LP;
3040 EV_TEMP[IEV_TAB_TEMP][4] = c * pc * zcv / ETOT_LP;
3043 EV_TEMP[IEV_TAB_TEMP][2] = pc * xcv;
3044 EV_TEMP[IEV_TAB_TEMP][3] = pc * ycv;
3045 EV_TEMP[IEV_TAB_TEMP][4] = pc * zcv;
3047 G4double VXOUT=0.,VYOUT=0.,VZOUT=0.;
3049 EV_TEMP[IEV_TAB_TEMP][2],EV_TEMP[IEV_TAB_TEMP][3],
3050 EV_TEMP[IEV_TAB_TEMP][4],
3051 &VXOUT,&VYOUT,&VZOUT);
3052 EV_TEMP[IEV_TAB_TEMP][2] = VXOUT;
3053 EV_TEMP[IEV_TAB_TEMP][3] = VYOUT;
3054 EV_TEMP[IEV_TAB_TEMP][4] = VZOUT;
3057 G4double v2 = std::pow(EV_TEMP[IEV_TAB_TEMP][2],2.) +
3058 std::pow(EV_TEMP[IEV_TAB_TEMP][3],2.) +
3059 std::pow(EV_TEMP[IEV_TAB_TEMP][4],2.);
3060 G4double gamma = 1.0/std::sqrt(1.0 - v2 / (c*c));
3061 G4double etot_lp = amoins*mu * gamma;
3062 pxeva = pxeva - EV_TEMP[IEV_TAB_TEMP][2] * etot_lp /
c;
3063 pyeva = pyeva - EV_TEMP[IEV_TAB_TEMP][3] * etot_lp /
c;
3064 pleva = pleva - EV_TEMP[IEV_TAB_TEMP][4] * etot_lp /
c;
3067 pxeva = pxeva - EV_TEMP[IEV_TAB_TEMP][2];
3068 pyeva = pyeva - EV_TEMP[IEV_TAB_TEMP][3];
3069 pleva = pleva - EV_TEMP[IEV_TAB_TEMP][4];
3071 G4double pteva = std::sqrt(pxeva*pxeva + pyeva*pyeva);
3073 G4double etot = std::sqrt ( pleva*pleva + pteva*pteva + af*af * mu2 );
3074 vxeva = c * pxeva / etot;
3075 vyeva = c * pyeva / etot;
3076 vleva = c * pleva / etot;
3077 IEV_TAB_TEMP = IEV_TAB_TEMP + 1;
3080 if(twon==1){
goto secondneutron;}
3083 if (zf < 3. || (ff == 1) || (fgamma == 1) || (fimf==1)) {
3095 (*tkeimf_par) = tkeimf;
3096 (*mtota_par) = mtota;
3097 (*vleva_par) = vleva;
3098 (*vxeva_par) = vxeva;
3099 (*vyeva_par) = vyeva;
3102 (*inttype_par) = inttype;
3103 (*iev_tab_temp_par)= IEV_TAB_TEMP;
3105 (*NbLam0_par) = NbLam0;
3109 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 *problamb0_par,
G4double *ptotl_par,
G4double *sn_par,
G4double *sbp_par,
G4double *sbd_par,
G4double *sbt_par,
G4double *sbhe_par,
G4double *sba_par,
G4double *slamb0_par,
G4double *ecn_par,
G4double *ecp_par,
G4double *ecd_par,
G4double *ect_par,
G4double *eche_par,
G4double *eca_par,
G4double *ecg_par,
G4double *eclamb0_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 *jprflamb0_par,
G4double *tsum_par,
G4int NbLam0)
3120 G4double problamb0 = (*problamb0_par);
3273 G4int choice_fisspart = 0;
3382 mglw(a-1.0,zprf,&ma1z);
3383 mglw(a-1.0,zprf-1.0,&ma1z1);
3384 mglw(a-2.0,zprf-1.0,&ma2z1);
3385 mglw(a-3.0,zprf-1.0,&ma3z1);
3386 mglw(a-3.0,zprf-2.0,&ma3z2);
3387 mglw(a-4.0,zprf-2.0,&ma4z2);
3390 if((a-1.)==3.0 && (zprf-1.0)==2.0) ma1z1=-7.7181660;
3391 if((a-1.)==4.0 && (zprf-1.0)==2.0) ma1z1=-28.295992;
3396 sd = ma2z1 - maz - 2.2246;
3397 st = ma3z1 - maz - 8.481977;
3398 she = ma3z2 - maz - 7.7181660;
3399 sa = ma4z2 - maz - 28.295992;
3429 if (zprf <= 1.0e0 || a <= 1.0e0 || (a-zprf) < 0.0) {
3439 if (zprf <= 1.0e0 || a <= 2.0e0 || (a-zprf) < 1.0) {
3449 if (zprf <= 1.0e0 || a <= 3.0e0 || (a-zprf) < 2.0) {
3459 if (a-4.0<=0.0 || zprf<=2.0 || (a-zprf)<2.0) {
3469 if (a-3.0 <= 0.0 || zprf<=2.0 || (a-zprf)<1.0) {
3483 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);
3494 barfit(k,k+j,il,&sbfis,&segs,&selmax);
3500 ef = ef*(4.5114-2.2687*(a-zprf)/zprf);
3502 ef = ef*(3.3931-1.5338*(a-zprf)/zprf);
3506 if((a-zprf)/zprf>1.52)ef=ef*(1.1222-0.10886*(a-zprf)/zprf)-0.1;
3508 if(k>=94&&k<=98&&j<158){
3510 if(
mod(j,2)==0&&
mod(k,2)==0){
3511 if(k>=94){ef = ef-(11.54108*(a-zprf)/zprf-18.074);}
3514 if(
mod(j,2)==1&&
mod(k,2)==1){
3515 if(k>=95){ef = ef-(14.567*(a-zprf)/zprf-23.266);}
3518 if(
mod(j,2)==0&&
mod(k,2)==1){
3519 if(j>=144){ef = ef-(13.662*(a-zprf)/zprf-21.656);}
3522 if(
mod(j,2)==1&&
mod(k,2)==0){
3523 if(j>=144){ef = ef-(13.662*(a-zprf)/zprf-21.656);}
3534 if (ef < 0.0)ef = 0.0;
3540 ef = ef + 0.51*(1115.-938.+sn-slamb0)/std::pow(a,2./3.);
3573 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);
3574 erot = jprf * jprf * 197.328 * 197.328 /(2. * iinert);
3577 bsbkbc(a,zprf,&bscn,&bkcn,&bccn);
3580 densniv(a,zprf,ee,0.0,&densg,bshell,bscn,bkcn,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprf,0,&qrcn);
3632 lorb(a,a-1.,jprf,ee-sn,&dlout,&sdlout);
3634 if(IDjprf==1) djprf = 0.0;
3635 jprfn = jprf + djprf;
3641 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);
3642 erotn = jprfn * jprfn * 197.328 * 197.328 /(2. * iinert);
3643 bsbkbc(a-1.,zprf,&bs,&bk,&bc);
3646 densniv(a-1.0,zprf,ee,sn,&densn,bshell, bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprfn,0,&qr);
3656 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ_NEUT CALLED MORE THAN 100 TIMES" << std::endl;
3665 if(ecn<=0.0)
goto dir1234;
3685 lorb(a,a-1.,jprf,ee-sbp,&dlout,&sdlout);
3687 if(IDjprf==1) djprf = 0.0;
3688 jprfp = jprf + djprf;
3694 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);
3695 erotp = jprfp * jprfp * 197.328 * 197.328 /(2. * iinert);
3697 bsbkbc(a-1.,zprf-1.,&bs,&bk,&bc);
3700 densniv(a-1.0,zprf-1.0,ee,sbp,&densp,bshell,bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprfp,0,&qr);
3710 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3719 if(ecp<=0.0)
goto dir1235;
3722 ecp = 2.0 * pt +
bp;
3736 if ((in >= 2) && (iz >= 2)) {
3740 lorb(a,a-2.,jprf,ee-sbd,&dlout,&sdlout);
3742 if(IDjprf==1) djprf = 0.0;
3743 jprfd = jprf + djprf;
3749 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);
3750 erotd = jprfd * jprfd * 197.328 * 197.328 /(2. * iinert);
3752 bsbkbc(a-2.,zprf-1.,&bs,&bk,&bc);
3755 densniv(a-2.0,zprf-1.0e0,ee,sbd,&densd,bshell,bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprfd,0,&qr);
3766 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3775 if(ecd<=0.0)
goto dir1236;
3778 ecd = 2.0 * dt + bd;
3792 if ((in >= 3) && (iz >= 2)) {
3796 lorb(a,a-3.,jprf,ee-sbt,&dlout,&sdlout);
3798 if(IDjprf==1) djprf = 0.0;
3799 jprft = jprf + djprf;
3805 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);
3806 erott = jprft * jprft * 197.328 * 197.328 /(2. * iinert);
3808 bsbkbc(a-3.,zprf-1.,&bs,&bk,&bc);
3811 densniv(a-3.0,zprf-1.0,ee,sbt,&denst,bshell,bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprft,0,&qr);
3822 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3831 if(ect<=0.0)
goto dir1237;
3834 ect = 2.0 * tt + bt;
3848 if ((in >= 3) && (iz >= 3)) {
3852 lorb(a,a-4.,jprf,ee-sba,&dlout,&sdlout);
3854 if(IDjprf==1) djprf = 0.0;
3855 jprfa = jprf + djprf;
3861 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);
3862 erota = jprfa * jprfa * 197.328 * 197.328 /(2. * iinert);
3864 bsbkbc(a-4.,zprf-2.,&bs,&bk,&bc);
3867 densniv(a-4.0,zprf-2.0,ee,sba,&densa,bshell,bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprfa,0,&qr);
3878 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3887 if(eca<=0.0)
goto dir1238;
3890 eca = 2.0 * at + ba;
3904 if ((in >= 2) && (iz >= 3)) {
3908 lorb(a,a-3.,jprf,ee-sbhe,&dlout,&sdlout);
3910 if(IDjprf==1) djprf = 0.0;
3911 jprfhe = jprf + djprf;
3917 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);
3918 erothe = jprfhe * jprfhe * 197.328 * 197.328 /(2. * iinert);
3920 bsbkbc(a-3.,zprf-2.,&bs,&bk,&bc);
3923 densniv(a-3.0,zprf-2.0,ee,sbhe,&denshe,bshell,bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprfhe,0,&qr);
3934 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3943 if(eche<=0.0)
goto dir1239;
3946 eche = 2.0 * het + bhe;
3962 if (in >= 2 && NbLam0>0) {
3966 lorb(a,a-1.,jprf,ee-slamb0,&dlout,&sdlout);
3968 if(IDjprf==1) djprf = 0.0;
3969 jprflamb0 = jprf + djprf;
3975 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);
3976 erotlamb0 = jprflamb0 * jprflamb0 * 197.328 * 197.328 /(2. * iinert);
3977 bsbkbc(a-1.,zprf,&bs,&bk,&bc);
3980 densniv(a-1.0,zprf,ee,slamb0,&denslamb0,bshell, bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprflamb0,0,&qr);
3990 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ_NEUT CALLED MORE THAN 100 TIMES" << std::endl;
3993 if(eclamb0>(ee-slamb0)){
3994 if((ee-slamb0)<rlamb0t)
3995 eclamb0 = ee-slamb0;
3999 if(eclamb0<=0.0)
goto dir1240;
4001 eclamb0 = 2.0 * lamb0t;
4029 gn =
width(a,zprf,1.0,0.0,nt,0.0,sn,ee-erotn)* densn/densg;
4034 gp =
width(a,zprf,1.0,1.0,pt,bp,sbp,ee-erotp)*densp/densg*
pen(a, 1.0, omegap, pt);
4039 gd =
width(a,zprf,2.0,1.0,dt,bd,sbd,ee-erotd)*densd/densg*
pen(a, 2.0, omegad, dt);
4044 gt =
width(a,zprf,3.0,1.0,tt,bt,sbt,ee-erott)*denst/densg*
pen(a, 3.0, omegat, tt);
4049 ghe =
width(a,zprf,3.0,2.0,het,bhe,sbhe,ee-erothe) * denshe/densg*
pen(a, 3.0, omegahe, het);
4054 ga =
width(a,zprf,4.0,2.0,at,ba,sba,ee-erota) * densa/densg*
pen(a, 4.0, omegaa, at);
4059 glamb0 =
width(a,zprf,1.0,-2.0,lamb0t,0.0,slamb0,ee-erotlamb0)* denslamb0/densg;
4067 G4int izcn=0,incn=0,inmin=0,inmax=0,inmi=0,inma=0;
4070 if(fimf_allowed==0 || zprf<=5.0 || a<=7.0){
4088 inmin =
max(inmin,incn-inma);
4089 inmax =
min(inmax,incn-inmi);
4091 inmax =
max(inmax,inmin);
4093 for(
G4int iaimf=izimf+inmin;iaimf<=izimf+inmax;iaimf++){
4095 if(aimf>=a || zimf>=zprf){
4107 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));
4109 erot = jprf * jprf * 197.328 * 197.328 /(2.0 * iinert);
4112 if(densg==0.0 || ee < (sbimf + erot)){
4118 densniv(a,zprf,ee,sbimf,&densimf,0.0,bsimf,1.0,&timf,0,0,defbetimf,&ecor,jprf,2,&qr);
4120 imfarg = (sbimf+erotcn-erot)/timf;
4121 if(imfarg > 200.0) imfarg = 200.0;
4131 width_imf =
width(a,zprf,aimf,zimf,timf,bimf,sbimf,ee-erot)*std::exp(-imfarg)*qr/qrcn;
4134 gimf3 = gimf3 + width_imf;
4148 inmin =
max(inmin,incn-inma);
4149 inmax =
min(inmax,incn-inmi);
4151 inmax =
max(inmax,inmin);
4153 for(
G4int iaimf=izimf+inmin;iaimf<=izimf+inmax;iaimf++){
4155 if(aimf>=a || zimf>=zprf){
4167 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));
4169 erot = jprf * jprf * 197.328 * 197.328 /(2.0 * iinert);
4172 if(densg==0.0 || ee < (sbimf + erot)){
4178 densniv(a,zprf,ee,sbimf,&densimf,0.0,bsimf,1.0,&timf,0,0,defbetimf,&ecor,jprf,2,&qr);
4180 imfarg = (sbimf+erotcn-erot)/timf;
4181 if(imfarg > 200.0) imfarg = 200.0;
4190 width_imf =
width(a,zprf,aimf,zimf,timf,bimf,sbimf,ee-erot)*std::exp(-imfarg)*qr/qrcn;
4193 gimf5 = gimf5 + width_imf;
4198 if(gimf3<=0.0 || gimf5<=0.0){
4204 b_imf = (std::log10(gimf3) - std::log10(gimf5))/(std::log10(3.0)-std::log10(5.0));
4206 if(b_imf >= -1.01) b_imf = -1.01;
4207 if(b_imf <= -100.0) {
4214 a_imf = gimf3 / std::pow(3.0,b_imf);
4215 gimf = a_imf * ( std::pow(zprf,b_imf+1.0) - std::pow(3.0,b_imf+1.0)) /(b_imf + 1.0);
4219 if(gimf < 1.
e-10) gimf = 0.0;
4226 pa = (
ald->
av)*a + (
ald->
as)*std::pow(a,2./3.) + (
ald->
ak)*std::pow(a,1./3.);
4227 gamma = 2.5 * pa * std::pow(a,-4./3.);
4233 gtemp = 17.60/(std::pow(a,0.699) * std::sqrt(gfactor));
4236 gg = 0.624e-9*std::pow(a,1.6)*std::pow(gtemp,5.);
4246 gsum = ga + ghe + gd + gt + gp + gn + gimf + gg + glamb0;
4259 if(
fiss->
ifis==0 || (zprf*zprf/a<=22.74 && zprf<60.)){
4267 fission_width(zprf,a,ee,bssp,bksp,ef,y,&gf,&temp,jprf,0,1,
fiss->
optcol,
fiss->
optshp,densg);
4287 gtotal = ga + ghe + gp + gd + gt + gn + gg +gimf + gf + glamb0;
4306 probhe = ghe/gtotal;
4309 probimf = gimf/gtotal;
4310 problamb0 = glamb0/gtotal;
4323 fomega_sp(a,y,&mfcd,&omegasp,&homegasp);
4327 fomega_gs(a,zprf,&k1,&omegags,&homegags);
4340 part_fiss(
fiss->
bet,gsum,gf,y,tauc,ts1,tsum, &choice_fisspart,zprf,a,ft,&t_lapse,&gff);
4344 tsum = tsum + t_lapse;
4347 if(choice_fisspart==2){
4363 gtotal=ga + ghe + gp + gd + gt + gn + gimf + gg + glamb0;
4382 probhe = ghe/gtotal;
4385 probimf = gimf/gtotal;
4386 problamb0 = glamb0/gtotal;
4392 gtotal = ga + ghe + gp + gd + gt + gn + gg + gimf + glamb0;
4410 probhe = ghe/gtotal;
4413 probimf = gimf/gtotal;
4414 problamb0 = glamb0/gtotal;
4418 ptotl = probp+probd+probt+probn+probhe+proba+probg+probimf+probf+problamb0;
4424 (*probp_par) = probp;
4425 (*probd_par) = probd;
4426 (*probt_par) = probt;
4427 (*probn_par) = probn;
4428 (*probhe_par) = probhe;
4429 (*proba_par) = proba;
4430 (*probg_par) = probg;
4431 (*probimf_par) = probimf;
4432 (*problamb0_par) = problamb0;
4433 (*probf_par) = probf;
4434 (*ptotl_par) = ptotl;
4441 (*slamb0_par) = slamb0;
4454 (*eclamb0_par) = eclamb0;
4462 (*jprfn_par) = jprfn;
4463 (*jprfp_par) = jprfp;
4464 (*jprfd_par) = jprfd;
4465 (*jprft_par) = jprft;
4466 (*jprfhe_par) = jprfhe;
4467 (*jprfa_par) = jprfa;
4468 (*jprflamb0_par) = jprflamb0;
4473 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)
4568 G4double pi6 = std::pow(3.1415926535,2) / 6.0;
4580 if(afp<=20) BSHELLCT = 0.0;
4613 pa = (
ald->
av)*a + (
ald->
as)*std::pow(a,(2.e0/3.e0)) + (
ald->
ak)*std::pow(a,(1.e0/3.e0));
4615 pa = (
ald->
av)*a + (
ald->
as)*bsin*std::pow(a,(2.e0/3.e0)) + (
ald->
ak)*bkin*std::pow(a,(1.e0/3.e0));
4617 gamma = 2.5 * pa * std::pow(a,-4.0/3.0);
4622 if(ifis==0&&bs!=1.0){
4625 if(ponq>700.0) ponq = 700.0;
4626 bs = 1.0/(1.0+std::exp(-ponq)) + 1.0/(1.0+std::exp(ponq)) * bsin;
4627 bk = 1.0/(1.0+std::exp(-ponq)) + 1.0/(1.0+std::exp(ponq)) * bkin;
4632 pa = (
ald->
av)*a + (
ald->
as)*std::pow(a,(2.e0/3.e0)) + (
ald->
ak)*std::pow(a,(1.e0/3.e0));
4635 pa = (
ald->
av)*a + (
ald->
as)*bs*std::pow(a,(2.e0/3.e0)) + (
ald->
ak)*bk*std::pow(a,(1.e0/3.e0));
4638 gamma = 2.5 * pa * std::pow(a,-4.0/3.0);
4644 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));
4664 deltpp = -0.25e0* std::pow((delta0/std::sqrt(a)),2) * pa /pi6 + 22.34e0*std::pow(a,-0.464)-0.235;
4668 e=e-(0.285+11.17*std::pow(a,-0.464)-0.390-0.00058*
a);
4672 e=e-(22.34*std::pow(a,-0.464)-0.235);
4696 ponfe = -2.5*pa*e*std::pow(a,(-4.0/3.0));
4698 if (ponfe < -700.0) {
4701 fe = 1.0 - std::exp(ponfe);
4704 he = 1.0 - std::pow((1.0 - e/ecr),2);
4711 fecor = e + deltau*fe + deltpp*he;
4718 y1 = std::sqrt(pa*fecor);
4719 for(
G4int j = 0; j < 5; j++) {
4720 y2 = pa*fecor*(1.e0-std::exp(-y1));
4725 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;
4728 y11 = std::sqrt(pa*ecor1);
4729 for(
G4int j = 0; j < 7; j++) {
4730 y21 = pa*ecor1*(1.0-std::exp(-y11));
4731 y11 = std::sqrt(y21);
4735 fdens = fdens*std::pow((y01/y0),1.5);
4736 ftemp = ftemp*std::pow((y01/y0),1.5);
4740 ponniv = 2.0*std::sqrt(pa*fecor);
4741 if (ponniv > 700.0) {
4745 fdens = 0.1477045 * std::exp(ponniv)/(std::pow(pa,0.25)*std::pow(fecor,1.25));
4746 ftemp = std::sqrt(fecor/pa);
4752 if(IOPTCT==0)
goto densniv100;
4753 tempct = 17.60/( std::pow(a,0.699) * std::sqrt(1.+gamma*BSHELLCT));
4764 if (IPARITE == 1) { e0 = 0.285+11.17*std::pow(a,-0.464) - 0.390-0.00058*
a;}
4766 if (IPARITE == 2) { e0 = 22.34*std::pow(a,-0.464)-0.235;}
4768 if (IPARITE == 0){ e0 = 0.0;}
4770 ponniv = (ein-e0)/tempct;
4771 if(ifis!=1) ponniv =
max(0.0,(ein-e0)/tempct);
4772 if(ponniv>700.0){ ponniv = 700.0;}
4773 densct = std::exp(ponniv)/tempct*std::exp(0.079*BSHELLCT/tempct);
4777 if(elim>=ecr&&densfm<=densct){
4785 if(elim>=ecr&&tfm>=tempct){
4793 ponniv = (ein)/tempct;
4794 if(ponniv>700.0){ ponniv = 700.0;}
4795 densct = std::exp(ponniv)/tempct;
4797 if(ein>=ecr && densfm<=densct){
4807 if(ein>=ecr && tfm>=tempct){
4822 ftemp = 17.60/( std::pow(a,0.699) * std::sqrt(1.0+gamma*BSHELLCT));
4835 fnorm = std::pow(1.16,2)*931.49*1.e-2/(9.0* std::pow(6.582122,2));
4837 if(ifis==0 || ifis==2){
4843 fp_per = 0.4*std::pow(a,5.0/3.0)*fnorm*(1.0+0.50*defbet*std::sqrt(5.0/(4.0*pi)));
4844 fp_par = 0.40*std::pow(a,5.0/3.0)*fnorm*(1.0-defbet*std::sqrt(5.0/(4.0*pi)));
4853 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));
4855 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));
4862 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);
4863 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);
4867 if(fp_par<0.0)fp_par=0.0;
4868 if(fp_per<0.0)fp_per=0.0;
4870 sig_per = std::sqrt(fp_per * ftemp);
4871 sig_par = std::sqrt(fp_par * ftemp);
4873 sigma2 = sig_per*sig_per + sig_par*sig_par;
4874 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));
4875 erot = jprf*jprf/(2.0*std::sqrt(fp_par*fp_par+fp_per*fp_per));
4879 qrot(z,a,defbet,sig_per,fecor-erot,&fqr);
4885 fdens = fdens * fqr *jfact;
4887 if(fdens<1e-300)fdens=0.0;
4924 G4int distn,distz,ndist, zdist;
4925 G4int nmn[8]= {2, 8, 14, 20, 28, 50, 82, 126};
4926 G4int nmz[8]= {2, 8, 14, 20, 28, 50, 82, 126};
4941 for(
G4int i =0;i<8;i++){
4942 ndist = std::fabs(
idnint(
n) - nmn[i]);
4943 if(ndist < distn) distn = ndist;
4944 zdist = std::fabs(
idnint(z) - nmz[i]);
4945 if(zdist < distz) distz = zdist;
4951 bet = 0.022 + 0.003*dn + 0.002*
dz;
4953 sig = 75.0*std::pow(bet,2.) * sig;
4957 ponq = (u - ucr)/dcr;
4965 (*qr) = 1.0/(1.0 + std::exp(ponq)) * (sig - 1.0) + 1.0;
4986 for(
G4int i = 2; i <
n; i++) {
5012 if(ia==0)
return eflmacResult;
5015 G4double z = 0.0,
n = 0.0,
a = 0.0, av = 0.0, as = 0.0;
5016 G4double a0 = 0.0,
c1 = 0.0, c4 = 0.0, b1 = 0.0, b3 = 0.0;
5018 G4double r0 = 0.0, kf = 0.0, ks = 0.0;
5019 G4double kv = 0.0, rp = 0.0, ay = 0.0, aden = 0.0, x0 = 0.0, y0 = 0.0;
5020 G4double esq = 0.0, ael = 0.0, i = 0.0, e0 = 0.0;
5070 if(flag==1){
goto eflmac311;}
5080 c1 = 3.0/5.0*esq/r0;
5081 c4 = 5.0/4.0*std::pow((3.0/(2.0*pi)),(2.0/3.0)) *
c1;
5082 kf = std::pow((9.0*pi*z/(4.0*
a)),(1.0/3.0))/r0;
5084 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));
5087 x0 = r0 * std::pow(
a,(1.0/3.0)) / ay;
5088 y0 = r0 * std::pow(
a,(1.0/3.0)) / aden;
5090 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);
5092 b3 = 1.0 - 5.0/std::pow(y0,2) * (1.0 - 15.0/(8.0*y0) + 21.0/(8.0 * std::pow(y0,3))
5093 - 3.0/4.0 * (1.0 + 9.0/(2.0*y0) + 7.0/std::pow(y0,2)
5094 + 7.0/(2.0 * std::pow(y0,3))) * std::exp(-2.0*y0));
5098 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
5099 +
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))
5100 + ff*std::pow(z,2)/
a -ca*(
n-
z) - ael * std::pow(z,(2.39e0));
5107 if (in==iz && (
mod(in,2) == 1) && (
mod(iz,2) == 1) && in>0.) {
5123 e0 = 0.285+11.17*std::pow(
a,-0.464) -0.390-0.00058*(
a);
5129 e0 = 22.34*std::pow(
a,-0.464)-0.235;
5141 return eflmacResult;
5166 (*del) = -12.0/std::sqrt(a);
5169 (*del) = 12.0/std::sqrt(a);
5217 if (bet/(std::sqrt(2.0)*10.0*(homega/6.582122)) <= 1.0) {
5218 tauResult = std::log(10.0*ef/t)/(bet*1.0e21);
5221 tauResult = std::log(10.0*ef/t)/ (2.0*std::pow((10.0*homega/6.582122),2))*(bet*1.0
e-21);
5233 G4double rel = bet/(20.0*homega/6.582122);
5234 G4double cramResult = std::sqrt(1.0 + std::pow(rel,2)) - rel;
5237 if (cramResult > 1.0) {
5259 const G4int bsbkSize = 54;
5261 G4double bk[bsbkSize] = {0.0, 1.00000,1.00087,1.00352,1.00799,1.01433,1.02265,1.03306,
5262 1.04576,1.06099,1.07910,1.10056,1.12603,1.15651,1.19348,
5263 1.23915,1.29590,1.35951,1.41013,1.44103,1.46026,1.47339,
5264 1.48308,1.49068,1.49692,1.50226,1.50694,1.51114,1.51502,
5265 1.51864,1.52208,1.52539,1.52861,1.53177,1.53490,1.53803,
5266 1.54117,1.54473,1.54762,1.55096,1.55440,1.55798,1.56173,
5267 1.56567,1.56980,1.57413,1.57860,1.58301,1.58688,1.58688,
5268 1.58688,1.58740,1.58740, 0.0};
5270 G4double bs[bsbkSize] = {0.0, 1.00000,1.00086,1.00338,1.00750,1.01319,
5271 1.02044,1.02927,1.03974,
5272 1.05195,1.06604,1.08224,1.10085,1.12229,1.14717,1.17623,1.20963,
5273 1.24296,1.26532,1.27619,1.28126,1.28362,1.28458,1.28477,1.28450,
5274 1.28394,1.28320,1.28235,1.28141,1.28042,1.27941,1.27837,1.27732,
5275 1.27627,1.27522,1.27418,1.27314,1.27210,1.27108,1.27006,1.26906,
5276 1.26806,1.26707,1.26610,1.26514,1.26418,1.26325,1.26233,1.26147,
5277 1.26147,1.26147,1.25992,1.25992, 0.0};
5279 i =
idint(y/(2.0
e-02)) + 1;
5281 if((i + 1) >= bsbkSize) {
5289 bipolResult = bs[i] + (bs[i+1] - bs[i])/2.0
e-02 * (y - 2.0
e-02*(i - 1));
5292 bipolResult = bk[i] + (bk[i+1] - bk[i])/2.0
e-02 * (y - 2.0
e-02*(i - 1));
5307 ES0 = 20.760*std::pow(AF,2.0/3.0);
5310 MR02 = std::pow(AF,5.0/3.0)*1.0340*0.010*1.175*1.175;
5312 (*MFCD) = MR02 * 3.0/10.0*(1.0+3.0*
Y);
5314 OMEGA = std::sqrt(ES0/MR02)*std::sqrt(8.0/3.0*Y*(1.0+304.0*Y/255.0));
5316 HOMEGA = 6.58122*OMEGA/10.0;
5331 G4double OMEGA,HOMEGA,MR02,MINERT,
C,fk1;
5333 MR02 = std::pow(AF,5.0/3.0)*1.0340*0.01*1.175*1.175;
5334 MINERT = 3.*MR02/10.0;
5335 C = 17.9439*(1.-1.7826*std::pow((AF-2.0*ZF)/AF,2));
5336 fk1 = 0.4*C*std::pow(AF,2.0/3.0)-0.1464*std::pow(ZF,2)/std::pow(AF,1./3.);
5337 OMEGA = std::sqrt(fk1/MINERT);
5338 HOMEGA = 6.58122*OMEGA/10.0;
5373 BARR = 1.345 * Z1 * Z2 / RMAX;
5375 OMEGA = 4.5 / 197.3287;
5459 for(
G4int init_i = 0; init_i < 7; init_i++) {
5463 for(
G4int init_i = 0; init_i < 10; init_i++) {
5467 G4double a = 0.0,
z = 0.0, amin = 0.0, amax = 0.0, amin2 = 0.0;
5468 G4double amax2 = 0.0, aa = 0.0,
zz = 0.0, bfis = 0.0;
5469 G4double bfis0 = 0.0, ell = 0.0, el = 0.0, egs = 0.0, el80 = 0.0, el20 = 0.0;
5470 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;
5471 G4double aj = 0.0, ak = 0.0, a1 = 0.0, a2 = 0.0;
5473 G4int i = 0, j = 0,
k = 0,
m = 0;
5476 G4double emncof[4][5] = {{-9.01100e+2,-1.40818e+3, 2.77000e+3,-7.06695e+2, 8.89867e+2},
5477 {1.35355e+4,-2.03847e+4, 1.09384e+4,-4.86297e+3,-6.18603e+2},
5478 {-3.26367e+3, 1.62447e+3, 1.36856e+3, 1.31731e+3, 1.53372e+2},
5479 {7.48863e+3,-1.21581e+4, 5.50281e+3,-1.33630e+3, 5.05367e-2}};
5481 G4double elmcof[4][5] = {{1.84542e+3,-5.64002e+3, 5.66730e+3,-3.15150e+3, 9.54160e+2},
5482 {-2.24577e+3, 8.56133e+3,-9.67348e+3, 5.81744e+3,-1.86997e+3},
5483 {2.79772e+3,-8.73073e+3, 9.19706e+3,-4.91900e+3, 1.37283e+3},
5484 {-3.01866e+1, 1.41161e+3,-2.85919e+3, 2.13016e+3,-6.49072e+2}};
5486 G4double emxcof[4][6] = {{9.43596e4,-2.241997e5,2.223237e5,-1.324408e5,4.68922e4,-8.83568e3},
5487 {-1.655827e5,4.062365e5,-4.236128e5,2.66837e5,-9.93242e4,1.90644e4},
5488 {1.705447e5,-4.032e5,3.970312e5,-2.313704e5,7.81147e4,-1.322775e4},
5489 {-9.274555e4,2.278093e5,-2.422225e5,1.55431e5,-5.78742e4,9.97505e3}};
5491 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},
5492 {-1.13269453e+6, 2.97764590e+6,-4.54326326e+6, 3.00464870e+6, -1.44989274e+6,-1.02026610e+5, 6.27959815e+4},
5493 {1.37543304e+6,-3.65808988e+6, 5.47798999e+6,-3.78109283e+6, 1.84131765e+6, 1.53669695e+4,-6.96817834e+4},
5494 {-8.56559835e+5, 2.48872266e+6,-4.07349128e+6, 3.12835899e+6, -1.62394090e+6, 1.19797378e+5, 4.25737058e+4},
5495 {3.28723311e+5,-1.09892175e+6, 2.03997269e+6,-1.77185718e+6, 9.96051545e+5,-1.53305699e+5,-1.12982954e+4},
5496 {4.15850238e+4, 7.29653408e+4,-4.93776346e+5, 6.01254680e+5, -4.01308292e+5, 9.65968391e+4,-3.49596027e+3},
5497 {-1.82751044e+5, 3.91386300e+5,-3.03639248e+5, 1.15782417e+5, -4.24399280e+3,-6.11477247e+3, 3.66982647e+2}};
5499 const G4int sizex = 5;
5500 const G4int sizey = 6;
5501 const G4int sizez = 4;
5503 G4double egscof[sizey][sizey][sizez];
5505 G4double egs1[sizey][sizex] = {{1.927813e5, 7.666859e5, 6.628436e5, 1.586504e5,-7.786476e3},
5506 {-4.499687e5,-1.784644e6,-1.546968e6,-4.020658e5,-3.929522e3},
5507 {4.667741e5, 1.849838e6, 1.641313e6, 5.229787e5, 5.928137e4},
5508 {-3.017927e5,-1.206483e6,-1.124685e6,-4.478641e5,-8.682323e4},
5509 {1.226517e5, 5.015667e5, 5.032605e5, 2.404477e5, 5.603301e4},
5510 {-1.752824e4,-7.411621e4,-7.989019e4,-4.175486e4,-1.024194e4}};
5512 G4double egs2[sizey][sizex] = {{-6.459162e5,-2.903581e6,-3.048551e6,-1.004411e6,-6.558220e4},
5513 {1.469853e6, 6.564615e6, 6.843078e6, 2.280839e6, 1.802023e5},
5514 {-1.435116e6,-6.322470e6,-6.531834e6,-2.298744e6,-2.639612e5},
5515 {8.665296e5, 3.769159e6, 3.899685e6, 1.520520e6, 2.498728e5},
5516 {-3.302885e5,-1.429313e6,-1.512075e6,-6.744828e5,-1.398771e5},
5517 {4.958167e4, 2.178202e5, 2.400617e5, 1.167815e5, 2.663901e4}};
5519 G4double egs3[sizey][sizex] = {{3.117030e5, 1.195474e6, 9.036289e5, 6.876190e4,-6.814556e4},
5520 {-7.394913e5,-2.826468e6,-2.152757e6,-2.459553e5, 1.101414e5},
5521 {7.918994e5, 3.030439e6, 2.412611e6, 5.228065e5, 8.542465e3},
5522 {-5.421004e5,-2.102672e6,-1.813959e6,-6.251700e5,-1.184348e5},
5523 {2.370771e5, 9.459043e5, 9.026235e5, 4.116799e5, 1.001348e5},
5524 {-4.227664e4,-1.738756e5,-1.795906e5,-9.292141e4,-2.397528e4}};
5526 G4double egs4[sizey][sizex] = {{-1.072763e5,-5.973532e5,-6.151814e5, 7.371898e4, 1.255490e5},
5527 {2.298769e5, 1.265001e6, 1.252798e6,-2.306276e5,-2.845824e5},
5528 {-2.093664e5,-1.100874e6,-1.009313e6, 2.705945e5, 2.506562e5},
5529 {1.274613e5, 6.190307e5, 5.262822e5,-1.336039e5,-1.115865e5},
5530 {-5.715764e4,-2.560989e5,-2.228781e5,-3.222789e3, 1.575670e4},
5531 {1.189447e4, 5.161815e4, 4.870290e4, 1.266808e4, 2.069603e3}};
5533 for(i = 0; i < sizey; i++) {
5534 for(j = 0; j < sizex; j++) {
5535 egscof[i][j][0] = egs1[i][j];
5536 egscof[i][j][1] = egs2[i][j];
5537 egscof[i][j][2] = egs3[i][j];
5538 egscof[i][j][3] = egs4[i][j];
5543 if (iz < 19 || iz > 111) {
5547 if(iz > 102 && il > 0) {
5554 amin= 1.2e0*
z + 0.01e0*
z*
z;
5555 amax= 5.8e0*z - 0.024e0*z*
z;
5557 if(a < amin || a > amax) {
5569 for(i = 0; i < 7; i++) {
5570 for(j = 0; j < 7; j++) {
5571 bfis0=bfis0+elzcof[j][i]*pz[i]*pa[j];
5583 amin2 = 1.4e0*z + 0.009e0*z*
z;
5584 amax2 = 20.e0 + 3.0e0*
z;
5586 if((a < amin2-5.e0 || a > amax2+10.e0) && il > 0) {
5596 for(i = 0; i < 4; i++) {
5597 for(j = 0; j < 5; j++) {
5598 el80 = el80 + elmcof[i][j]*pz[j]*pa[i];
5599 el20 = el20 + emncof[i][j]*pz[j]*pa[i];
5610 for(i = 0; i < 4; i++) {
5611 for(j = 0; j < 6; j++) {
5612 elmax = elmax + emxcof[i][j]*pz[j]*pa[i];
5623 x = sel20/(*selmax);
5624 y = sel80/(*selmax);
5628 q = 0.2/(std::pow(sel20,2)*std::pow(sel80,2)*(sel20-sel80));
5629 qa = q*(4.0*std::pow(sel80,3) - std::pow(sel20,3));
5630 qb = -q*(4.0*std::pow(sel80,2) - std::pow(sel20,2));
5631 bfis = bfis*(1.0 + qa*std::pow(el,2) + qb*std::pow(el,3));
5635 aj = (-20.0*std::pow(
x,5) + 25.e0*std::pow(
x,4) - 4.0)*std::pow((
y-1.0),2)*
y*
y;
5636 ak = (-20.0*std::pow(y,5) + 25.0*std::pow(y,4) - 1.0) * std::pow((
x-1.0),2)*
x*
x;
5637 q = 0.2/(std::pow((y-x)*((1.0-x)*(1.0-y)*x*y),2));
5638 qa = q*(aj*y - ak*
x);
5639 qb = -q*(aj*(2.0*y + 1.0) - ak*(2.0*x + 1.0));
5641 a1 = 4.0*std::pow(z,5) - 5.0*std::pow(z,4) + 1.0;
5642 a2 = qa*(2.e0*z + 1.e0);
5643 bfis=bfis*(a1 + (z - 1.e0)*(a2 + qb*z)*z*z*(z - 1.e0));
5650 if(el > (*selmax)) {
5656 if(el > (*selmax)) {
5660 for(
k = 0;
k < 4;
k++) {
5661 for(l = 0; l < 6; l++) {
5662 for(
m = 0;
m < 5;
m++) {
5663 egs = egs + egscof[l][
m][
k]*pz[l]*pa[
k]*pl[2*
m];
5709 ferf=-
gammp(0.5,x*x);
5711 ferf=
gammp(0.5,x*x);;
5721 if(x<0.0 || a<=0.0)std::cout <<
"G4Abla::gammp = bad arguments in gammp" << std::endl;
5723 gser(&gamser,a,x,gln);
5726 gcf(&gammcf,a,x,gln);
5745 for(
G4int i=1;i<=itmax;i++){
5749 if(std::fabs(d)<fpmin)d=fpmin;
5751 if(std::fabs(c)<fpmin)c=fpmin;
5755 if(std::fabs(del-1.)<
eps)
goto dir1;
5757 std::cout <<
"a too large, ITMAX too small in gcf" << std::endl;
5759 fgammcf=std::exp(-x+a*std::log(x)-gln)*
h;
5772 if(x<0.)std::cout <<
"G4Abla::gser = x < 0 in gser" << std::endl;
5783 if(std::fabs(del)<std::fabs(sum)*eps)
goto dir1;
5785 std::cout <<
"a too large, ITMAX too small in gser" << std::endl;
5787 fgamser=sum*std::exp(-x+a*std::log(x)-gln);
5795 G4double cof[6]={76.18009172947146,-86.50532032941677,24.01409824083091,
5796 -1.231739572450155,0.1208650973866179e-2,-0.5395239384953e-5};
5802 tmp=(x+0.5)*std::log(tmp)-
tmp;
5803 ser=1.000000000190015;
5804 for(
G4int j=0;j<6;j++){
5809 return fgammln=tmp+std::log(stp*ser/x);
5817 return (E*std::exp(-E));
5823 return (1.0 - (E + 1.0) * std::exp(-E));
5839 const G4int pSize = 101;
5857 for(i = 1; i <= 99; i++) {
5861 if (std::fabs(
f(x) -
G4double(i)/100.0) < 1
e-5) {
5886 x = (p[i] - p[i-1])*(y*100 - i) + p[i];
5903 if(ii <= 0 || jj < 0) {
5914 fpace2=fpace2/1000.;
5916 if(
pace->
dm[ii][jj] == 0.) {
5921 guet(&a, &z, &fpace2);
5922 fpace2=fpace2-ii*931.5;
5923 fpace2=fpace2/1000.;
5942 const G4int qrows = 50;
5943 const G4int qcols = 70;
5945 for(
G4int init_i = 0; init_i < qrows; init_i++) {
5946 for(
G4int init_j = 0; init_j < qcols; init_j++) {
5947 q[init_i][init_j] = 0.0;
5988 G4double aux=1.+(9.*xjj/4./qq/x13);
5990 G4double ee4=avol*xx+asur*(std::pow(xx,.666))+ac*x13+azer;
5991 G4double tota = ee1 + ee2 + ee3 + ee4;
5992 find = 939.55*xneu+938.77*zz - tota;
6004 const G4double fmp = 938.27231,fmn=939.56563,fml=1115.683;
6010 for(
G4int i=0;i<IMULTBU;i++){
6026 G4double gamma = std::sqrt(1.0 - v2 / (c*c));
6027 G4double avvmass = iz*fmp + (ia-iz-is)*fmn + is*fml +
eflmac(ia,iz,0,3);
6037 for(
G4int i=0;i<IEV_TAB;i++){
6053 G4double gamma = std::sqrt(1.0 - v2 / (c*c));
6054 G4double avvmass = iz*fmp + (ia-iz-is)*fmn + is*fml +
eflmac(ia,iz,0,3);
6067 G4double gamma = std::sqrt(1.0 - v2 / (c*c));
6162 fractpart = std::modf(number, &intpart);
6167 if(fractpart < 0.5) {
6168 return G4int(std::floor(number));
6171 return G4int(std::ceil(number));
6175 if(fractpart < -0.5) {
6176 return G4int(std::floor(number));
6179 return G4int(std::ceil(number));
6183 return G4int(std::floor(number));
6192 mylocaltime = localtime(&mytime);
6195 return(mylocaltime->tm_hour*60*60 + mylocaltime->tm_min*60 + mylocaltime->tm_sec);
6223 if(x-std::floor(x) <= std::ceil(x)-x)
6234 if(x-std::floor(x) <= std::ceil(x)-x)
6235 value =
G4int(std::floor(x));
6237 value =
G4int(std::ceil(x));
6244 if(x-std::floor(x) <= std::ceil(x)-x)
6245 return G4int(std::floor(x));
6247 return G4int(std::ceil(x));
6252 if(a < b && a < c) {
6255 if(b < a && b < c) {
6258 if(c < a && c < b) {
6278 G4int IZPART,IAPART,NMOTHER;
6281 G4double INT1,INT2,INT3,AKONST,EARG,R0,MPARTNER;
6284 G4double PAR_A1=0.,PAR_B1=0.,FACT=1.;
6299 NMOTHER =
idnint(AMOTHER-ZMOTHER);
6308 HBAR = 6.582122e-22;
6313 APARTNER = AMOTHER - APART;
6314 MPARTNER = APARTNER * 931.49 /
C2;
6317 if(IAPART==1&&IZPART==0){
6319 MPART = 939.56 /
C2;
6320 if(idlamb0==1)MPART = 1115.683 /
C2;
6322 if(IAPART==1&&IZPART==1){
6324 MPART = 938.27 /
C2;
6327 if(IAPART==2&&IZPART==0){
6329 MPART = 2.*939.56 /
C2;
6331 if(IAPART==2&&IZPART==1){
6333 MPART = 1876.10 /
C2;
6335 if(IAPART==3&&IZPART==1){
6337 MPART = 2809.39 /
C2;
6339 if(IAPART==3&&IZPART==2){
6341 MPART = 2809.37 /
C2;
6343 if(IAPART==4&&IZPART==2){
6345 MPART = 3728.35 /
C2;
6349 MPART = APART * 931.49 /
C2;
6359 MU = MPARTNER * MPART / (MPARTNER + MPART);
6363 RGEOM = R0 * (std::pow(APART,1.0/3.0)+std::pow(AMOTHER-APART,1.0/3.0));
6366 AKONST = HBAR*std::sqrt(1.0 / MU);
6369 BKONST = MPART / ( PI * PI * HBAR * HBAR);
6373 INT1 = 2.0 * std::pow(TEMP,3.) / (2.0 * TEMP +
B);
6375 ARG = std::sqrt(B/TEMP);
6376 EARG = (
erf(ARG) - 1.0);
6377 if(
std::abs(EARG)<1.e-9) EARG = 0.0;
6379 INT2 = 0.5 * std::sqrt(PI) * std::pow(TEMP,3.0/2.0);
6382 if(AEXP>700.0) AEXP = 700.0;
6383 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;
6384 if(INT2<0.0) INT2 = 0.0;
6387 if(EARG==0.0) INT2 = 0.0;
6390 INT3 = 2.0*TEMP*TEMP*TEMP / (2.0*TEMP*TEMP + 4.0*B*TEMP + B*
B);
6392 if(IZPART<-1.0&&ZMOTHER<151.0){
6396 fwidth = PI * BKONST * G * std::sqrt((RGEOM * RGEOM * INT1 + 2.0 * AKONST * RGEOM * INT2 + AKONST * AKONST * INT3) * RGEOM * RGEOM * INT1);
6399 fwidth = PI * BKONST * G *(RGEOM * RGEOM * INT1 + 2.0 * AKONST * RGEOM * INT2 + AKONST * AKONST * INT3);
6407 PAR_A1=std::exp(2.302585*0.2083*std::exp(-0.01548472*AMOTHER))-0.05;
6408 PAR_B1 = 0.59939389 + 0.00915657 * AMOTHER;
6410 if(AMOTHER>154.0&&AMOTHER<195.0){
6411 PAR_A1=1.0086961-8.629e-5*AMOTHER;
6412 PAR_B1 = 1.5329331 + 0.00302074 * AMOTHER;
6414 if(AMOTHER>194.0&&AMOTHER<208.0){
6415 PAR_A1=9.8356347-0.09294663*AMOTHER+2.441e-4*AMOTHER*AMOTHER;
6416 PAR_B1 = 7.7701987 - 0.02897401 * AMOTHER;
6418 if(AMOTHER>207.0&&AMOTHER<228.0){
6419 PAR_A1=15.107385-0.12414415*AMOTHER+2.7222e-4*AMOTHER*AMOTHER;
6420 PAR_B1=-64.078009+0.56813179*AMOTHER-0.00121078*AMOTHER*AMOTHER;
6423 if(
mod(NMOTHER,2)==0&&NMOTHER>147.){
6424 PAR_A1 = 2.0*(0.9389118 + 6.4559e-5 * AMOTHER);
6426 if(
mod(NMOTHER,2)==1)PAR_A1 = 3.0*(0.9389118 + 6.4559
e-5 * AMOTHER);
6428 PAR_B1 = 2.1507177 + 0.00146119 * AMOTHER;
6434 FACT = std::exp((2.302585*PAR_A1*std::exp(-PAR_B1*(EXC-SB))));
6435 if(FACT<1.0) FACT = 1.0;
6436 if(IZPART<-1.&&ZMOTHER<151.0){
6438 fwidth = fwidth / std::sqrt(FACT);
6440 fwidth = fwidth / FACT;
6445 std::cout <<
"LOOK IN PARTICLE_WIDTH!" << std::endl;
6446 std::cout <<
"ACN,APART :"<< AMOTHER << APART << std::endl;
6447 std::cout <<
"EXC,TEMP,B,SB :" << EXC <<
" " << TEMP <<
" " << B <<
" " << SB << std::endl;
6448 std::cout <<
"INTi, i=1-3 :" << INT1 <<
" " << INT2 <<
" " << INT3 << std::endl;
6449 std::cout <<
" " << std::endl;
6466 MU = (A -
ap) * ap / A;
6470 HO = 197.3287 * omega;
6475 fpen=std::pow(10.0,4.
e-4*std::pow(T/(HO*HO*std::pow(MU,0.25)),-4.3/2.3026));
6493 (*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;
6495 (*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;
6539 G4double DEFO_INIT,OMEGA,HOMEGA,OMEGA_GS,HOMEGA_GS,K1,MFCD;
6540 G4double BET1,XACT,SIGMA_SQR,W_EXP,XB,NORM,SIGMA_SQR_INF,W_INFIN,W;
6541 G4double FUNC_TRANS,LOG_SLOPE_INF,LOG_SLOPE_ABS;
6548 fomega_gs(AF,ZF,&K1,&OMEGA_GS,&HOMEGA_GS);
6552 if((bet*bet)>4.0*OMEGA_GS*OMEGA_GS){
6553 BET1=std::sqrt(bet*bet-4.0*OMEGA_GS*OMEGA_GS);
6558 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)));
6561 XACT = DEFO_INIT *std::exp(-0.5*(bet-BET1)*1.e21*(TIME-T_0));
6566 BET1=std::sqrt(4.0*OMEGA_GS*OMEGA_GS-bet*bet);
6567 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));
6568 XACT = DEFO_INIT*std::cos(0.5*BET1*1.e21*(TIME-T_0))*std::exp(-bet*1.e21*(TIME-T_0));
6574 XB = 7./3.*Y-938./765.*Y*Y+9.499768*Y*Y*Y-8.050944*Y*Y*Y*
Y;
6579 NORM = 1./std::sqrt(2.*PI*SIGMA_SQR);
6581 W_EXP = -1.*(XB - XACT)*(XB - XACT)/(2.0 * SIGMA_SQR);
6582 if(W_EXP<(-708.0) ) W_EXP = -708.0;
6583 W = NORM * std::exp( W_EXP ) * FT / (K1 * SIGMA_SQR);
6590 SIGMA_SQR_INF = FT/K1;
6591 W_EXP = -XB*XB/(2.0 * SIGMA_SQR_INF);
6592 if(W_EXP<(-708.0))W_EXP = -708.0;
6593 W_INFIN = std::exp(W_EXP)/std::sqrt(2.0*PI*SIGMA_SQR_INF);
6594 FUNC_TRANS = W / W_INFIN;
6599 LOG_SLOPE_INF =
cram(bet,HOMEGA)*bet*MFCD*OMEGA/FT;
6600 LOG_SLOPE_ABS = (XB-XACT)/SIGMA_SQR-XB/SIGMA_SQR_INF+
cram(bet,HOMEGA)*bet*MFCD*OMEGA/FT;
6602 FUNC_TRANS = FUNC_TRANS * LOG_SLOPE_ABS/LOG_SLOPE_INF;
6608 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)
6667 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;
6685 if(BET*BET>4.0*OMEGA*OMEGA){
6692 t_0 = BET*1.e21*HBAR*HBAR/(4.*HOMEGA*FT)/16.;
6695 if(((2.*FT-HOMEGA/16.)>0.000001) && BET>0.0){
6700 t_0 = (std::log(2.*FT/(2.*FT-HOMEGA/16.)))/(BET*1.e21);
6712 STEP_LENGTH = 1.5*TAUF/50.;
6717 BEGIN_TIME = TSUM + t_0;
6719 if(BEGIN_TIME<0.0) std::cout <<
"CURRENT TIME < 0" << BEGIN_TIME << std::endl;
6721 if(BEGIN_TIME<1.50*TAUF){
6722 LOC_TIME_BEGIN = BEGIN_TIME;
6724 while((LOC_TIME_BEGIN<1.5*TAUF)&&fchoice==0){
6726 LOC_TIME_END = LOC_TIME_BEGIN + STEP_LENGTH;
6729 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;
6731 fGF_LOC = fGF_LOC * GF;
6741 LAMBDA = 1.0/TS1 + 1.0/TS2;
6747 REAC_PROB = std::exp(-1.0*STEP_LENGTH*LAMBDA);
6752 FISS_PROB = fGF_LOC / (fGF_LOC+GP);
6763 LOC_TIME_BEGIN = LOC_TIME_END;
6766 fT_LAPSE = LOC_TIME_END - BEGIN_TIME;
6773 FISS_PROB = GF / (GF+GP);
6783 LAMBDA = 1./TS1 + 1./TS2;
6805 (*T_LAPSE)=fT_LAPSE;
6818 G4double MFCD,OMEGA,HOMEGA1,HOMEGA2=0.,GFTUN;
6819 G4double E1,E2,EXP_FACT,CORR_FUNCT,FACT1,FACT2,FACT3;
6827 if(
mod(IN,2)==0&&
mod(IZ,2)==0){
6830 EE = EE - 12.0/std::sqrt(A);
6834 if(
mod(IN,2)==1&&
mod(IZ,2)==1){
6838 if(
mod(IN,2)==1&&
mod(IZ,2)==0){
6842 if(
mod(IN,2)==0&&
mod(IZ,2)==1){
6846 E1 = EF + HOMEGA1/2.0/PI*std::log(HOMEGA1*(2.0*PI+HOMEGA2)/4.0/PI/PI);
6848 E2 = EF + HOMEGA2/(2.0*PI)*std::log(1.0+2.0*PI/HOMEGA2);
6855 EXP_FACT = (EE-EF)/(HOMEGA2/(2.0*PI));
6856 if(EXP_FACT>700.0) EXP_FACT = 700.0;
6857 CORR_FUNCT = HOMEGA1 * (1.0-1.0/(1.0+std::exp(EXP_FACT)));
6858 if(
mod(IN,2)==0&&
mod(IZ,2)==0){
6859 CORR_FUNCT = HOMEGA1 * (1.0-1.0/(1.0+std::exp(EXP_FACT)));
6862 FACT1 = HOMEGA1/(2.0*PI*TEMP+HOMEGA1);
6863 FACT2 = (2.0*PI/(2.0*PI+HOMEGA2)-HOMEGA1*(2.0*PI+HOMEGA2)/4.0/PI/PI)/(E2-E1);
6864 FACT3 = HOMEGA2/(2.0*PI*TEMP-HOMEGA2);
6867 GFTUN = FACT1*(std::exp(EE/TEMP)*std::exp(2.0*PI*(EE-EF)/HOMEGA1)-std::exp(-2.0*PI*EF/HOMEGA1));
6870 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));
6872 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));
6875 GFTUN = GFTUN/std::exp(EE/TEMP)*DENSF*ENH_FACT/DENSG/2.0/PI;
6876 GFTUN = GFTUN * CORR_FUNCT;
6881 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)
6884 G4double FNORM,MASS_ASYM_SADD_B,FP_PER,FP_PAR,SIG_PER_SP,SIG_PAR_SP;
6885 G4double Z2OVERA,ftemp,fgf,DENSF,ECOR,EROT,qr;
6886 G4double DCR,UCR,ENH_FACTA,ENH_FACTB,ENH_FACT,PONFE;
6891 Z2OVERA = ZPRF * ZPRF /
A;
6894 if((ZPRF<=55.0) || (FF_ALLOWED==0)){
6904 densniv(A,ZPRF,EE,EF,&DENSF,0.0,BS,BK,&ftemp,OPTSHP,0,Y,&ECOR,JPR,1,&qr);
6907 fgf= DENSF/DENSG/PI/2.0*ftemp;
6919 FNORM = 1.2*1.2 * 931.49 * 1.e-2 / (9.0 * 6.582122*6.582122);
6922 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));
6927 if(Z2OVERA<=30.0) FP_PER = 6.50;
6930 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));
6931 if(FP_PAR<0.0) FP_PAR = 0.0;
6933 EROT = JPR * JPR / (2.0 * std::sqrt(FP_PAR*FP_PAR + FP_PER*FP_PER));
6934 if(IEROT==1) EROT = 0.0;
6937 SIG_PER_SP = std::sqrt(FP_PER * ftemp);
6939 if(SIG_PER_SP<1.0) SIG_PER_SP = 1.0;
6942 SIG_PAR_SP = std::sqrt(FP_PAR * ftemp);
6946 MASS_ASYM_SADD_B = 2.0;
6948 MASS_ASYM_SADD_B = 1.0;
6952 if(Z2OVERA>35.&&Z2OVERA<=(110.*110./298.0)){
6954 ENH_FACTA = std::sqrt(8.0*PI) * SIG_PER_SP*SIG_PER_SP * SIG_PAR_SP;
6956 ENH_FACTB = MASS_ASYM_SADD_B * SIG_PER_SP*SIG_PER_SP;
6958 ENH_FACT = ENH_FACTA * ENH_FACTB / (ENH_FACTA + ENH_FACTB);
6962 ENH_FACT = MASS_ASYM_SADD_B*SIG_PER_SP*SIG_PER_SP;
6965 ENH_FACT = std::sqrt(8.0*PI) * SIG_PER_SP*SIG_PER_SP* SIG_PAR_SP;
6970 PONFE = (ECOR-UCR-EROT)/DCR;
6971 if(PONFE>700.) PONFE = 700.0;
6973 ENH_FACT = 1.0/(1.0+std::exp(PONFE))*ENH_FACT+1.0;
6975 if(ENH_FACT<1.0)ENH_FACT = 1.0;
6976 fgf= DENSF/DENSG/PI/2.0*ftemp*ENH_FACT;
6980 fgf=
tunnelling(A,ZPRF,Y,EE,EF,ftemp,DENSG,DENSF,ENH_FACT);
6992 G4double AFRAGMENT,S4FINAL,ALEVDENS;
6993 G4double THETA_MOTHER,THETA_ORBITAL;
7008 if (EEFINAL<=0.01) EEFINAL = 0.01;
7009 AFRAGMENT = AMOTHER - ADAUGHTER;
7010 ALEVDENS = 0.073*AMOTHER + 0.095*std::pow(AMOTHER,2.0/3.0);
7011 S4FINAL = ALEVDENS * EEFINAL;
7012 if(S4FINAL <= 0.0 || S4FINAL > 100000.){
7013 std::cout<<
"S4FINAL:" << S4FINAL << ALEVDENS << EEFINAL <<
idnint(AMOTHER) <<
idnint(AFRAGMENT) << std::endl;
7015 THETA_MOTHER = 0.0111 * std::pow(AMOTHER,1.66667);
7016 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);
7018 *LORBITAL = -1.* THETA_ORBITAL * (LMOTHER / THETA_MOTHER + std::sqrt(S4FINAL) /(ALEVDENS*LMOTHER));
7020 *SIGMA_LORBITAL = std::sqrt(std::sqrt(S4FINAL) * THETA_ORBITAL / ALEVDENS);
7048 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;
7049 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;
7050 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;
7051 G4double PI = 3.141592653589793238;
7053 G4double SDWprevious=0,SUMDW_TOT=0,SUM_Z=0,
X=0,SUMDW_N_TOT=0,XX=0;
7061 IZIMFMAX =
idnint(ZCN / 2.0);
7064 std::cout <<
"CHARGE_IMF line 46" << std::endl;
7065 std::cout <<
"Problem: IZIMFMAX < 3 " << std::endl;
7066 std::cout <<
"ZCN,IZIMFMAX," << ZCN <<
"," << IZIMFMAX << std::endl;
7074 bsbkbc(ACN,ZCN,&BS,&BK,&BC);
7076 densniv(ACN,ZCN,EE,0.0,&DENSCN,BSHELL,BS,BK,&TEMPCN,0,0,DEFBET,&ECOR,JPRF,0,&QRCN);
7078 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);
7079 EROTCN = JPRF * JPRF * 197.328 * 197.328 /(2. * IINERT);
7081 for(IZIMF=3;IZIMF<=IZIMFMAX;IZIMF++){
7090 INIMFMI =
max(1,INIMFMI-2);
7093 INCN =
idnint(ACN) - IZCN;
7097 INMI =
max(1,INMI-2);
7098 INMIN =
max(INIMFMI,INCN-INMA);
7099 INNMAX =
min(INIMFMA,INCN-INMI);
7101 ILIMMAX =
max(INNMAX,INMIN);
7104 for(
G4int INIMF=INMIN;INIMF<=ILIMMAX;INIMF++){
7105 IAIMF = IZIMF + INIMF;
7106 DW[IZIMF][IAIMF] = 0.0;
7107 AIMF_1 = 1.0*(IAIMF);
7110 mglms(ACN-AIMF_1,ZCN-ZIMF_1,OPTSHPIMF,&MARES);
7111 mglms(AIMF_1,ZIMF_1,OPTSHPIMF,&MAIMF);
7112 mglms(ACN,ZCN,OPTSHPIMF,&MAZ);
7116 SSBIMF[IZIMF][IAIMF] = 1.e37;
7119 SSBIMF[IZIMF][IAIMF] = MAIMF + MARES - MAZ + fBIMF;
7120 BBIMF[IZIMF][IAIMF] = fBIMF;
7126 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));
7128 EROT = JPRF * JPRF * 197.328 * 197.328 /(2.0 * IINERT);
7131 if (EE<(SSBIMF[IZIMF][IAIMF]+EROT) || DENSCN<=0.0){
7140 densniv(ACN,ZCN,EE,SSBIMF[IZIMF][IAIMF],&DENSIMF,0.0,BSIMF,1.0,&fTIMF,0,0,DEFBETIMF,&ECOR,JPRF,2,&QR);
7141 IMFARG = (SSBIMF[IZIMF][IAIMF]+EROTCN-EROT)/fTIMF;
7142 if(IMFARG>200.0) IMFARG = 200.0;
7144 WIDTH1 =
width(ACN,ZCN,AIMF_1,ZIMF_1,fTIMF,fBIMF,SSBIMF[IZIMF][IAIMF],EE-EROT);
7146 WIDTH_IMF = WIDTH1 * std::exp(-IMFARG) * QR / QRCN;
7149 std::cout <<
"GAMMA_IMF=0 -> LOOK IN GAMMA_IMF CALCULATIONS!" << std::endl;
7150 std::cout <<
"ACN,ZCN,AIMF,ZIMF:" <<
idnint(ACN) <<
"," <<
idnint(ZCN) <<
"," <<
idnint(AIMF_1) <<
"," <<
idnint(ZIMF_1) << std::endl;
7151 std::cout <<
"SSBIMF,TIMF :" << SSBIMF[IZIMF][IAIMF] <<
"," << fTIMF << std::endl;
7152 std::cout <<
"DEXP(-IMFARG) = " << std::exp(-IMFARG) << std::endl;
7153 std::cout <<
"WIDTH1 =" << WIDTH1 << std::endl;
7157 SDW[IZIMF] = SDW[IZIMF] + WIDTH_IMF;
7159 DW[IZIMF][IAIMF] = WIDTH_IMF;
7167 SDWprevious = 1.e20;
7170 for(
G4int III_ZIMF=3;III_ZIMF<=IZIMFMAX;III_ZIMF++){
7172 if(SDW[III_ZIMF]==0.0){
7173 IZSTOP = III_ZIMF - 1;
7177 if(SDW[III_ZIMF]>SDWprevious){
7178 IZSTOP = III_ZIMF - 1;
7181 SDWprevious = SDW[III_ZIMF];
7193 A1PAR = std::log10(SDW[IZSTOP]/SDW[IZSTOP-2])/std::log10((1.0*IZSTOP)/(1.0*IZSTOP-2.0));
7194 A2PAR = std::log10(SDW[IZSTOP]) - A1PAR * std::log10(1.0*(IZSTOP));
7195 if(A2PAR>0.)A2PAR=-1.*A2PAR;
7196 if(A1PAR>0.)A1PAR=-1.*A1PAR;
7200 for(
G4int II_ZIMF = IZSTOP;II_ZIMF<=IZIMFMAX;II_ZIMF++){
7201 SDW[II_ZIMF] = std::pow(10.0,A2PAR) * std::pow(1.0*II_ZIMF,A1PAR);
7202 if(SDW[II_ZIMF]<0.0) SDW[II_ZIMF] = 0.0;
7209 for(
G4int I_ZIMF = 3;I_ZIMF<=IZIMFMAX;I_ZIMF++){
7210 SUMDW_TOT = SUMDW_TOT + SDW[I_ZIMF];
7213 std::cout <<
"*********************" << std::endl;
7214 std::cout <<
"IMF function" << std::endl;
7215 std::cout <<
"SUM of decay widths = " << SUMDW_TOT <<
" IZIMFMAX = " << IZIMFMAX << std::endl;
7216 std::cout <<
"IZSTOP = " << IZSTOP << std::endl;
7224 X =
haz(1)*SUMDW_TOT;
7231 for(
G4int IZ = 3;IZ<=IZIMFMAX;IZ++){
7232 SUM_Z = SUM_Z + SDW[IZ];
7245 INMINMEM =
max(1,INMINMEM-2);
7248 INMI =
max(1,INMI-2);
7251 INMINMEM =
max(INMINMEM,INCN-INMA);
7252 INMAXMEM =
min(INMAXMEM,INCN-INMI);
7254 INMAXMEM =
max(INMINMEM,INMAXMEM);
7258 for(
G4int IIINIMF = INMINMEM;IIINIMF<=INMAXMEM;IIINIMF++){
7259 IA = IZMEM + IIINIMF;
7260 if(IZMEM>=3&&IZMEM<=95&&IA>=4&&IA<=250){
7261 SUMDW_N_TOT = SUMDW_N_TOT + DW[IZMEM][IA];
7263 std::cout <<
"CHARGE IMF OUT OF RANGE" << IZMEM <<
", " << IA <<
", " <<
idnint(ACN) <<
", " <<
idnint(ZCN) <<
", " << TEMP << std::endl;
7267 XX =
haz(1)*SUMDW_N_TOT;
7270 for(
G4int IINIMF = INMINMEM;IINIMF<=INMAXMEM; IINIMF++){
7271 IIA = IZMEM + IINIMF;
7273 SUM_A = SUM_A + DW[IZMEM][IIA];
7282 NIMF = fAIMF - fZIMF;
7284 if((ACN-ZCN-NIMF)<=0.0 || (ZCN-fZIMF) <= 0.0){
7285 std::cout <<
"IMF Partner unstable:" << std::endl;
7286 std::cout <<
"System: Acn,Zcn,NCN:" << std::endl;
7287 std::cout <<
idnint(ACN) <<
", " <<
idnint(ZCN) <<
", " <<
idnint(ACN-ZCN) << std::endl;
7288 std::cout <<
"IMF: A,Z,N:" << std::endl;
7289 std::cout <<
idnint(fAIMF) <<
", " <<
idnint(fZIMF) <<
", " <<
idnint(fAIMF-fZIMF) << std::endl;
7290 std::cout <<
"Partner: A,Z,N:" << std::endl;
7291 std::cout <<
idnint(ACN-fAIMF) <<
", " <<
idnint(ZCN-fZIMF) <<
", " <<
idnint(ACN-ZCN-NIMF) << std::endl;
7292 std::cout <<
"----nmin,nmax" << INMINMEM <<
", " << INMAXMEM << std::endl;
7293 std::cout <<
"----- warning: Zimf=" << fZIMF <<
" Aimf=" << fAIMF << std::endl;
7294 std::cout <<
"----- look in subroutine IMF" << std::endl;
7295 std::cout <<
"ACN,ZCN,ZIMF,AIMF,temp,EE,JPRF::" << ACN <<
", " << ZCN <<
", " << fZIMF <<
", " << fAIMF <<
", " << TEMP <<
", " << EE <<
", " << JPRF << std::endl;
7296 std::cout <<
"-IZSTOP,IZIMFMAX:" << IZSTOP <<
", " << IZIMFMAX << std::endl;
7297 std::cout <<
"----X,SUM_Z,SUMDW_TOT:" <<
X <<
", " << SUM_Z <<
", " << SUMDW_TOT << std::endl;
7303 if(fZIMF>=ZCN || fAIMF>=ACN || fZIMF<=2 || fAIMF<=3){
7304 std::cout <<
"----nmin,nmax" << INMINMEM <<
", " << INMAXMEM << std::endl;
7305 std::cout <<
"----- warning: Zimf=" << fZIMF <<
" Aimf=" << fAIMF << std::endl;
7306 std::cout <<
"----- look in subroutine IMF" << std::endl;
7307 std::cout <<
"ACN,ZCN,ZIMF,AIMF,temp,EE,JPRF:" << ACN <<
", " << ZCN <<
", " << fZIMF <<
", " << fAIMF <<
", " << TEMP <<
", " << EE <<
", " << JPRF << std::endl;
7308 std::cout <<
"-IZSTOP,IZIMFMAX:" << IZSTOP <<
", " << IZIMFMAX << std::endl;
7309 std::cout <<
"----X,SUM_Z,SUMDW_TOT:" <<
X <<
", " << SUM_Z <<
", " << SUMDW_TOT << std::endl;
7310 for(
int III_ZIMF=3;III_ZIMF<=IZIMFMAX;III_ZIMF++)
7311 std::cout <<
"-**Z,SDW:" << III_ZIMF <<
", " << SDW[III_ZIMF] << std::endl;
7321 if((ZCN-fZIMF)<=0.0)std::cout <<
"CHARGE_IMF ZIMF > ZCN" << std::endl;
7322 if((ACN-fAIMF)<=0.0)std::cout <<
"CHARGE_IMF AIMF > ACN" << std::endl;
7327 EEDAUG = (EE - fSBIMF) * (ACN - fAIMF) / ACN;
7328 bsbkbc(ACN - fAIMF,ZCN-fZIMF,&BS,&BK,&BC);
7329 densniv(ACN-fAIMF,ZCN-fZIMF,EEDAUG,0.0,&DENSIMF,BSHELL,BS,BK,&fTIMF,0,0,DEFBET,&ECOR,0.0,0,&QR);
7332 std::cout <<
"----- warning: EE=" << EE <<
"," <<
" S+Bimf=" << fSBIMF << std::endl;
7333 std::cout <<
"----- look in subroutine IMF" << std::endl;
7334 std::cout <<
"IMF will be resampled" << std::endl;
7348 G4int VISOSTAB[191][2]={
7459 *nmin = VISOSTAB[z-1][0];
7460 *nmax = VISOSTAB[z-1][1];
7476 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, problamb0=0.0, ptotl = 0.0, tcn = 0.0;
7477 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, slamb0 = 0.0;
7478 G4double ecn = 0.0, ecp = 0.0, ecd = 0.0, ect = 0.0,eche = 0.0,eca = 0.0, ecg = 0.0, eclamb0 = 0.0,
bp = 0.0, bd = 0.0, bt = 0.0, bhe = 0.0, ba = 0.0;
7480 G4double xcv=0.,ycv=0.,zcv=0.,VXOUT=0.,VYOUT=0.,VZOUT=0.;
7482 G4double jprfn=0.0, jprfp=0.0, jprfd=0.0, jprft=0.0, jprfhe=0.0, jprfa=0.0, jprflamb0=0.0;
7483 G4double ctet1 = 0.0, stet1 = 0.0, phi1 = 0.0;
7493 G4int inttype=0,inum=0,gammadecay = 0, flamb0decay = 0;
7499 G4int NbLam0= (*NbLam0_par);
7503 const G4double mu2 = 931.494*931.494;
7523 a0 = 0.66482503 - 3.4678935 * std::exp(-0.0104002*ee);
7524 a1 = 5.6846e-04 + 0.00574515 * std::exp(-0.01114307*ee);
7525 tauf = (a0 + a1 * zf*zf/std::pow(af,0.3333333)) * tau0;
7528 direct(zf,af,ee,0.,&probp,&probd,&probt,&probn,&probhe,&proba,&probg,&probimf,&probf,&problamb0,&ptotl,
7529 &sn,&sbp,&sbd,&sbt,&sbhe,&sba,&slamb0,
7530 &ecn,&ecp,&ecd,&ect,&eche,&eca,&ecg,&eclamb0,
7531 &
bp,&bd,&bt,&bhe,&ba,&
sp,&sd,&st,&she,&sa,&ef,&ts1,inttype,inum,itest,&sortie,&tcn,
7532 &jprfn, &jprfp, &jprfd, &jprft, &jprfhe, &jprfa, &jprflamb0, &tsum, NbLam0);
7536 if(ptotl<=0.)
goto post100;
7540 if(emin>1e30)std::cout <<
"ERROR AT THE EXIT OF EVAPORA,E>1.D30,AF" << std::endl;
7547 pc = std::sqrt(std::pow((1.0 + ecn/9.3956
e2),2.) - 1.0) * 9.3956e2;
7551 else if(probp!=0.0){
7555 pc = std::sqrt(std::pow((1.0 + ecp/9.3827
e2),2.) - 1.0) * 9.3827e2;
7559 else if(probd!=0.0){
7563 pc = std::sqrt(std::pow((1.0 + ecd/1.875358
e3),2) - 1.0) * 1.875358e3;
7567 else if(probt!=0.0){
7571 pc = std::sqrt(std::pow((1.0 + ect/2.80828
e3),2) - 1.0) * 2.80828e3;
7575 else if(probhe!=0.0){
7578 epsiln = she + eche;
7579 pc = std::sqrt(std::pow((1.0 + eche/2.80826
e3),2) - 1.0) * 2.80826e3;
7583 else{
if(proba!=0.0){
7587 pc = std::sqrt(std::pow((1.0 + eca/3.72834
e3),2) - 1.0) * 3.72834e3;
7610 pc = std::sqrt(std::pow((1.0 + eca/3.72834
e3),2) - 1.0) * 3.72834e3;
7614 else if (
x < proba+probhe) {
7618 epsiln = she + eche;
7619 pc = std::sqrt(std::pow((1.0 + eche/2.80826
e3),2) - 1.0) * 2.80826e3;
7623 else if (
x < proba+probhe+probt) {
7628 pc = std::sqrt(std::pow((1.0 + ect/2.80828
e3),2) - 1.0) * 2.80828e3;
7632 else if (
x < proba+probhe+probt+probd) {
7637 pc = std::sqrt(std::pow((1.0 + ecd/1.875358
e3),2) - 1.0) * 1.875358e3;
7641 else if (
x < proba+probhe+probt+probd+probp) {
7646 pc = std::sqrt(std::pow((1.0 + ecp/9.3827
e2),2) - 1.0) * 9.3827e2;
7650 else if (
x < proba+probhe+probt+probd+probp+probn) {
7655 pc = std::sqrt(std::pow((1.0 + ecn/9.3956
e2),2.) - 1.0) * 9.3956e2;
7659 else if (
x < proba+probhe+probt+probd+probp+probn+problamb0) {
7663 epsiln = slamb0 + eclamb0;
7664 pc = std::sqrt(std::pow((1.0 + (eclamb0)/11.1568
e2),2.) - 1.0) * 11.1568e2;
7670 else if (
x < proba+probhe+probt+probd+probp+probn+problamb0+probg) {
7678 if(probp==0.0 && probn==0.0 && probd==0.0 && probt==0.0 && proba==0.0 && probhe==0.0 && problamb0==0.0 && probimf==0.0 && probf==0.0){
7689 if(gammadecay==1 && ee<=0.01+epsiln){
7698 if(ee<=0.01) ee = 0.010;
7700 if(af<2.5)
goto post100;
7716 ctet1 = 2.0*rnd - 1.0;
7717 stet1 = std::sqrt(1.0 - std::pow(ctet1,2));
7719 phi1 = rnd*2.0*3.141592654;
7720 xcv = stet1*std::cos(phi1);
7721 ycv = stet1*std::sin(phi1);
7726 G4double ETOT_LP = std::sqrt(pc*pc + amoins*amoins * mu2);
7727 if(flamb0decay==1)ETOT_LP = std::sqrt(pc*pc + 1115.683*1115.683);
7740 &VXOUT,&VYOUT,&VZOUT);
7750 G4double gamma = 1.0/std::sqrt(1.0 - v2 / (c*c));
7751 G4double etot_lp = amoins*mu * gamma;
7761 pteva = std::sqrt(pxeva*pxeva + pyeva*pyeva);
7763 etot = std::sqrt ( pleva*pleva + pteva*pteva + af*af * mu2 );
7764 vx_eva = c * pxeva / etot;
7765 vy_eva = c * pyeva / etot;
7766 vz_eva = c * pleva / etot;
7768 IEV_TAB_SSC = IEV_TAB_SSC +1;
7770 if(time<tauf)
goto post10;
7776 *E_scission_post = ee;
7777 *NbLam0_par = NbLam0;
7782 if(A<1.)
return (1.*H)/A*(10.68*A-21.27*std::pow(A,2./3.))*10.;
7783 return (1.*H)/A*(10.68*A-21.27*std::pow(A,2./3.));
7788 if(A<1.)
return 1.e38;
7792 if(Z==1 && A==4)
return 2.04;
7793 else if(Z==2 && A==4)
return 2.39;
7794 else if(Z==2 && A==5)
return 3.12;
7795 else if(Z==2 && A==6)
return 4.18;
7796 else if(Z==2 && A==7)
return 5.23;
7797 else if(Z==2 && A==8)
return 7.16;
7798 else if(Z==3 && A==6)
return 4.50;
7799 else if(Z==3 && A==7)
return 5.58;
7800 else if(Z==3 && A==8)
return 6.80;
7801 else if(Z==3 && A==9)
return 8.50;
7802 else if(Z==4 && A==7)
return 5.16;
7803 else if(Z==4 && A==8)
return 6.84;
7804 else if(Z==4 && A==9)
return 6.71;
7805 else if(Z==4 && A==10)
return 9.11;
7806 else if(Z==5 && A==9)
return 8.29;
7807 else if(Z==5 && A==10)
return 8.89;
7808 else if(Z==5 && A==11)
return 10.24;
7809 else if(Z==5 && A==12)
return 11.37;
7810 else if(Z==6 && A==12)
return 10.76;
7811 else if(Z==6 && A==13)
return 11.69;
7812 else if(Z==6 && A==14)
return 12.17;
7813 else if(Z==14 && A==28)
return 16.0;
7814 else if(Z==39 && A==89)
return 22.1;
7815 else if(Z==57 && A==139)
return 23.8;
7816 else if(Z==82 && A==208)
return 26.5;
7827 if(A<2 || Z<2)
return 0.;
7837 if(
mod(N,2) == 1 &&
mod(Z,2) == 1)
D = -12./std::sqrt(A);
7838 if(
mod(N,2) == 0 &&
mod(Z,2) == 0)
D = 12./std::sqrt(A);
7840 G4double deltanew = (1.-std::exp(-1.*A/
c))*
D;
7842 be= av*A-as*std::pow(A,2./3.)-ac*Z*(Z-1.)/std::pow(A,1./3.)-asym*(N-
Z)*(N-Z)/((1.+std::exp(-1.*A/
k))*A)+deltanew + ny*(0.0335*
my-26.7-48.7/std::pow(A,2.0/3.0));
7846 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)
7855 e =
dmin1(SBHE,SN,e);
7856 e =
dmin1(SBHE,SBA,e);
7877 *ECP = (-1.0)*SP + BP;
7894 *ECD = (-1.0)*SD + BD;
7911 *ECT = (-1.0)*ST + BT;
7928 *ECHE= (-1.0)*SHE + BHE;
7946 *ECA = (-1.0)*SA + BA;
8028 Delta_U1_shell_max = -2.45;
8034 Delta_U2_shell = -2.45;
8044 G4double Fwidth_asymm1,Fwidth_asymm2,Fwidth_symm;
8046 Fwidth_asymm1 = 0.65;
8047 Fwidth_asymm2 = 0.65;
8089 Friction_factor = 1.0;
8092 G4double cN_asymm1_shell, cN_asymm2_shell;
8093 G4double gamma,gamma_heavy1,gamma_heavy2;
8109 G4double Sasymm1=0.,Sasymm2=0.,Ssymm=0.,Ysum=0.,Yasymm=0.;
8111 G4double wNasymm1_saddle, wNasymm2_saddle, wNsymm_saddle;
8112 G4double wNasymm2_scission, wNsymm_scission;
8113 G4double wNasymm1, wNasymm2, wNsymm;
8136 G4double A_levdens_heavy1,A_levdens_heavy2;
8140 G4double epsilon_1_saddle,epsilon0_1_saddle;
8141 G4double epsilon_2_saddle,epsilon0_2_saddle,epsilon_symm_saddle;
8146 G4double E_eff1_saddle,E_eff2_saddle;
8147 G4double Epot0_mode1_saddle,Epot0_mode2_saddle,Epot0_symm_saddle;
8148 G4double Epot_mode1_saddle,Epot_mode2_saddle,Epot_symm_saddle;
8149 G4double E_defo,E_defo1,E_defo2,E_scission_pre,E_scission_post;
8155 G4double MassCurv_scis, MassCurv_sadd;
8157 G4double Nheavy1_shell,Nheavy2_shell;
8162 G4double Z_scission,N_scission,A_scission;
8166 G4double DSN132,Delta_U1_shell,E_eff0_saddle;
8167 G4int NbLam0= (*NbLam0_par);
8172 cN_asymm1_shell = 0.700 * N/
Z;
8173 cN_asymm2_shell = 0.040 * N/
Z;
8177 DSN132 = Nheavy1_in - N/Z * Zheavy1_in;
8178 Aheavy1 = Nheavy1_in + Zheavy1_in + 0.340 * DSN132;
8187 Delta_U1_shell = Delta_U1_shell_max + U1NZ_SLOPE *
std::abs(DSN132);
8188 Delta_U1_shell =
min(0.,Delta_U1_shell);
8192 Nheavy1 = N/A * Aheavy1;
8193 Aheavy2 = Nheavy2 * A/
N;
8197 A_levdens = A / xLevdens;
8198 gamma = A_levdens / (0.40 * std::pow(A,1.3333)) * FGAMMA;
8199 A_levdens_heavy1 = Aheavy1 / xLevdens;
8200 gamma_heavy1 = A_levdens_heavy1 / (0.40 * std::pow(Aheavy1,1.3333)) * FGAMMA * FGAMMA1;
8201 A_levdens_heavy2 = Aheavy2 / xLevdens;
8202 gamma_heavy2 = A_levdens_heavy2 / (0.40 * std::pow(Aheavy2,1.3333)) * FGAMMA;
8206 E_saddle_scission = (-24. + 0.02227 * Z*Z/std::pow(A,0.33333))*Friction_factor;
8207 E_saddle_scission =
max( 0.0, E_saddle_scission );
8213 Z2_over_A_eff = Z*Z/
A;
8215 if( Z2_over_A_eff< 34.0 )
8216 MassCurv_scis = std::pow(10., -1.093364 + 0.082933 * Z2_over_A_eff - 0.0002602 * Z2_over_A_eff*Z2_over_A_eff);
8218 MassCurv_scis = std::pow(10., 3.053536 - 0.056477 * Z2_over_A_eff+ 0.0002454 * Z2_over_A_eff*Z2_over_A_eff );
8223 MassCurv_sadd = X_s2s * MassCurv_scis;
8225 cN_symm = 8.0 / std::pow(N,2.) * MassCurv_scis;
8226 cN_symm_sadd = 8.0 / std::pow(N,2.) * MassCurv_sadd;
8228 Nheavy1_shell = Nheavy1;
8231 Nheavy1_eff = (cN_symm_sadd*Nsymm + cN_asymm1_shell *
8232 Uwash(E/A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1) *
8236 Uwash(E/A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1));
8238 Nheavy1_eff = (cN_symm_sadd*Nsymm +
8239 cN_asymm1_shell*Nheavy1_shell)
8244 Nheavy2_NZ = Nheavy2;
8245 Nheavy2_shell = Nheavy2_NZ;
8247 Nheavy2_eff = (cN_symm_sadd*Nsymm +
8249 Uwash(E/A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2) *
8253 Uwash(E/A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2));
8255 Nheavy2_eff = (cN_symm_sadd*Nsymm +
8256 cN_asymm2_shell*Nheavy2_shell)
8260 Delta_U1 = Delta_U1_shell + (Nheavy1_shell - Nheavy1_eff)*(Nheavy1_shell - Nheavy1_eff) * cN_asymm1_shell;
8261 Delta_U1 =
min(Delta_U1,0.0);
8262 Delta_U2 = Delta_U2_shell + (Nheavy2_shell - Nheavy2_eff)*(Nheavy2_shell - Nheavy2_eff) * cN_asymm2_shell;
8263 Delta_U2 =
min(Delta_U2,0.0);
8267 Epot0_mode1_saddle = (Nheavy1_eff-Nsymm)*(Nheavy1_eff-Nsymm) * cN_symm_sadd;
8268 Epot0_mode2_saddle = (Nheavy2_eff-Nsymm)*(Nheavy2_eff-Nsymm) * cN_symm_sadd;
8269 Epot0_symm_saddle = 0.0;
8273 Epot_mode1_saddle = Epot0_mode1_saddle + Delta_U1;
8274 Epot_mode2_saddle = Epot0_mode2_saddle + Delta_U2;
8275 Epot_symm_saddle = Epot0_symm_saddle;
8278 dUeff =
min( Epot_mode1_saddle, Epot_mode2_saddle);
8279 dUeff =
min( dUeff, Epot_symm_saddle);
8280 dUeff = dUeff - Epot_symm_saddle;
8289 epsilon0_1_saddle = Eld - Epot0_mode1_saddle;
8290 epsilon0_2_saddle = Eld - Epot0_mode2_saddle;
8293 epsilon_1_saddle = Eld - Epot_mode1_saddle;
8294 epsilon_2_saddle = Eld - Epot_mode2_saddle;
8296 epsilon_symm_saddle = Eld - Epot_symm_saddle;
8299 eexc1_saddle = epsilon_1_saddle;
8300 eexc2_saddle = epsilon_2_saddle;
8303 EEXC_MAX =
max( eexc1_saddle, eexc2_saddle);
8304 EEXC_MAX =
max( EEXC_MAX, Eld);
8307 epsilon_1_scission = Eld + E_saddle_scission - Epot_mode1_saddle;
8308 epsilon_2_scission = Eld + E_saddle_scission - Epot_mode2_saddle;
8311 epsilon_symm_scission = Eld + E_saddle_scission - Epot_symm_saddle;
8314 E_eff1_saddle = epsilon0_1_saddle - Delta_U1 *
8315 Uwash(epsilon_1_saddle/A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1);
8317 if( E_eff1_saddle < A_levdens * hbom1*hbom1)
8318 E_eff1_saddle = A_levdens * hbom1*hbom1;
8321 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_eff1_saddle) /
8323 Uwash(epsilon_1_saddle/A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)+
8326 E_eff2_saddle = epsilon0_2_saddle -
8328 Uwash(epsilon_2_saddle/A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2);
8330 if(E_eff2_saddle < A_levdens * hbom2*hbom2)
8331 E_eff2_saddle = A_levdens * hbom2*hbom2;
8334 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_eff2_saddle) /
8336 Uwash(epsilon_2_saddle/A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)+
8339 E_eff0_saddle = epsilon_symm_saddle;
8340 if(E_eff0_saddle < A_levdens * hbom3*hbom3)
8341 E_eff0_saddle = A_levdens * hbom3*hbom3;
8344 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_eff0_saddle) /
8347 if(epsilon_symm_scission > 0.0 ){
8348 E_HELP =
max(E_saddle_scission,epsilon_symm_scission);
8350 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*(E_HELP)) /
8354 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_saddle_scission) /
8361 if( E_saddle_scission == 0.0 ){
8362 wNasymm1_scission = wNasymm1_saddle;
8363 wNasymm2_scission = wNasymm2_saddle;
8365 if( Nheavy1_eff > 75.0 ){
8366 wNasymm1_scission = std::sqrt(21.0)*N/
A;
8367 wNasymm2_scission =
max( 12.8 - 1.0 *(92.0 - Nheavy2_eff),1.0)*N/
A;
8370 wNasymm1_scission = wNasymm1_saddle;
8371 wNasymm2_scission = wNasymm2_saddle;
8375 wNasymm1_scission =
max( wNasymm1_scission, wNasymm1_saddle );
8376 wNasymm2_scission =
max( wNasymm2_scission, wNasymm2_saddle );
8378 wNasymm1 = wNasymm1_scission * Fwidth_asymm1;
8379 wNasymm2 = wNasymm2_scission * Fwidth_asymm2;
8380 wNsymm = wNsymm_scission * Fwidth_symm;
8383 Aheavy1_eff = Nheavy1_eff * A/
N;
8384 Aheavy2_eff = Nheavy2_eff * A/
N;
8386 A_levdens_heavy1 = Aheavy1_eff / xLevdens;
8387 A_levdens_heavy2 = Aheavy2_eff / xLevdens;
8388 gamma_heavy1 = A_levdens_heavy1 / (0.40 * std::pow(Aheavy1_eff,1.3333)) * FGAMMA * FGAMMA1;
8389 gamma_heavy2 = A_levdens_heavy2 / (0.40 * std::pow(Aheavy2_eff,1.3333)) * FGAMMA;
8391 if( epsilon_symm_saddle < A_levdens * hbom3*hbom3)
8392 Ssymm = 2.0 * std::sqrt(A_levdens*A_levdens * hbom3*hbom3) +
8393 (epsilon_symm_saddle - A_levdens * hbom3*hbom3)/hbom3;
8395 Ssymm = 2.0 * std::sqrt(A_levdens*epsilon_symm_saddle);
8399 if( epsilon0_1_saddle < A_levdens * hbom1*hbom1 )
8400 Ssymm_mode1 = 2.0 * std::sqrt(A_levdens*A_levdens * hbom1*hbom1) +
8401 (epsilon0_1_saddle - A_levdens * hbom1*hbom1)/hbom1;
8403 Ssymm_mode1 = 2.0 * std::sqrt( A_levdens*epsilon0_1_saddle );
8405 if( epsilon0_2_saddle < A_levdens * hbom2*hbom2 )
8406 Ssymm_mode2 = 2.0 * std::sqrt(A_levdens*A_levdens * hbom2*hbom2) +
8407 (epsilon0_2_saddle - A_levdens * hbom2*hbom2)/hbom2;
8409 Ssymm_mode2 = 2.0 * std::sqrt(A_levdens*epsilon0_2_saddle);
8412 if( epsilon0_1_saddle -
8414 Uwash(epsilon_1_saddle/A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)
8415 < A_levdens * hbom1*hbom1 )
8416 Sasymm1 = 2.0 * std::sqrt( A_levdens*A_levdens * hbom1*hbom1 ) +
8417 (epsilon0_1_saddle - Delta_U1 *
8418 Uwash(epsilon_1_saddle/A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)
8419 - A_levdens * hbom1*hbom1)/hbom1;
8421 Sasymm1 = 2.0 *std::sqrt( A_levdens*(epsilon0_1_saddle - Delta_U1 *
8422 Uwash(epsilon_1_saddle/A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)));
8424 if( epsilon0_2_saddle -
8426 Uwash(epsilon_2_saddle/A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)
8427 < A_levdens * hbom2*hbom2 )
8428 Sasymm2 = 2.0 * std::sqrt( A_levdens*A_levdens * hbom2*hbom2 ) +
8429 (epsilon0_1_saddle-Delta_U1 *
8430 Uwash(epsilon_2_saddle/A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)
8431 - A_levdens * hbom2*hbom2)/hbom2;
8434 std::sqrt( A_levdens*(epsilon0_2_saddle - Delta_U2 *
8435 Uwash(epsilon_2_saddle/A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)));
8437 Yasymm1 = ( std::exp(Sasymm1 - Ssymm) - std::exp(Ssymm_mode1 - Ssymm) ) *
8438 wNasymm1_saddle / wNsymm_saddle * 2.0;
8440 Yasymm2 = ( std::exp(Sasymm2 - Ssymm) - std::exp(Ssymm_mode2 - Ssymm) ) *
8441 wNasymm2_saddle / wNsymm_saddle * 2.0;
8443 Ysum = Ysymm + Yasymm1 + Yasymm2;
8446 Ysymm = Ysymm / Ysum;
8447 Yasymm1 = Yasymm1 / Ysum;
8448 Yasymm2 = Yasymm2 / Ysum;
8449 Yasymm = Yasymm1 + Yasymm2;
8455 if( (epsilon_symm_saddle < epsilon_1_saddle) &&
8456 (epsilon_symm_saddle < epsilon_2_saddle) )
8459 if( epsilon_1_saddle < epsilon_2_saddle )
8467 r_e_o = std::pow(10.0,-0.0170 * (E_saddle_scission + Eld)*(E_saddle_scission + Eld));
8485 if( rmode < Yasymm1 )
8488 if( (rmode > Yasymm1) && (rmode < Yasymm) )
8497 N1mean = Nheavy1_eff;
8501 N1mean = Nheavy2_eff;
8517 while( N1r < 5.0 || N2r < 5.0 ){
8536 E_scission_pre =
max( epsilon_1_scission, 1.0 );
8538 if( N1mean > N*0.50 ){
8548 E_scission_pre =
max( epsilon_2_scission, 1.0 );
8549 if( N1mean > N*0.50 ){
8550 beta1 = (N1r - 92.0) * 0.030 + 0.60;
8555 beta1 =
max(beta1,beta1gs);
8556 beta2 = 1.0 - beta1;
8557 beta2 =
max(beta2,beta2gs);
8563 beta2 = (N2r -92.0) * 0.030 + 0.60;
8564 beta2 =
max(beta2,beta2gs);
8565 beta1 = 1.0 - beta2;
8566 beta1 =
max(beta1,beta1gs);
8581 beta =
max(0.177963+0.0153241*Zsymm-1.62037
e-4*Zsymm*Zsymm,betags);
8582 beta1 =
max(0.177963+0.0153241*Z1UCD-1.62037
e-4*Z1UCD*Z1UCD,beta1gs);
8583 beta2 =
max(0.177963+0.0153241*Z2UCD-1.62037
e-4*Z2UCD*Z2UCD,beta2gs);
8585 E_asym =
frldm( Z1UCD, N1r, beta1 ) +
8586 frldm( Z2UCD, N2r, beta2 ) +
8587 ecoul( Z1UCD, N1r, beta1, Z2UCD, N2r, beta2, 2.0 ) -
8588 2.0 *
frldm( Zsymm, Nsymm, beta ) -
8589 ecoul( Zsymm, Nsymm, beta, Zsymm, Nsymm, beta, 2.0 );
8590 E_scission_pre =
max( epsilon_symm_scission - E_asym, 1. );
8599 if(E_scission_pre>5. && NbLam0<1){
8601 &A_scission,&Z_scission,vx_eva_sc,vy_eva_sc,vz_eva_sc,&NbLam0);
8602 N_scission = A_scission - Z_scission;
8606 E_scission_post = E_scission_pre;
8607 N_scission = A_scission - Z_scission;
8613 N1r = N1r * N_scission /
N;
8614 N2r = N2r * N_scission /
N;
8615 Z1UCD = Z1UCD * Z_scission /
Z;
8616 Z2UCD = Z2UCD * Z_scission /
Z;
8654 CZ = (
frldm( Z1UCD-1.0, N1r+1.0, beta1 ) +
8655 frldm( Z2UCD+1.0, N2r-1.0, beta2 ) +
8656 frldm( Z1UCD+1.0, N1r-1.0, beta1 ) +
8657 frldm( Z2UCD-1.0, N2r+1.0, beta2 ) +
8658 ecoul( Z1UCD-1.0, N1r+1.0, beta1,
8659 Z2UCD+1.0, N2r-1.0, beta2, 2.0) +
8660 ecoul( Z1UCD+1.0, N1r-1.0, beta1,
8661 Z2UCD-1.0, N2r+1.0, beta2, 2.0) -
8662 2.0*
ecoul( Z1UCD, N1r, beta1, Z2UCD, N2r, beta2, 2.0) -
8663 2.0*
frldm( Z1UCD, N1r, beta1 ) -
8664 2.0*
frldm( Z2UCD, N2r, beta2) ) * 0.50;
8666 if(1.0/A_levdens*E_scission_post < 0.0)
8667 std::cout <<
"DSQRT 1 < 0" << A_levdens <<
" " << E_scission_post << std::endl;
8669 if(0.50 * std::sqrt(1.0/A_levdens*E_scission_post) / CZ < 0.0){
8670 std::cout <<
"DSQRT 2 < 0 " << CZ << std::endl;
8671 std::cout <<
"This event was not considered" << std::endl;
8675 ZA1width = std::sqrt(0.5*std::sqrt(1.0/A_levdens*E_scission_post)/CZ);
8686 ZA1width =
max(ZA1width,sigZmin);
8688 if(imode == 1 && cpol1 != 0.0){
8692 Z1rr = Z1UCD - cpol1 * A_scission/N_scission;
8698 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED" << std::endl;
8701 if ((
utilabs(Z1rr - Z1r) > 3.0*ZA1width) || Z1r<1.0)
goto fiss2801;
8704 if( imode == 2 && cpol2 != 0.0 ){
8708 Z1rr = Z1UCD - cpol2 * A_scission/N_scission;
8713 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED" << std::endl;
8716 if( (
utilabs(Z1rr - Z1r) > 3.0*ZA1width) || Z1r < 1.0 )
goto fiss2802;
8724 re1 =
frldm( Z1UCD-1.0, N1r+1.0, beta1 ) +
8725 frldm( Z2UCD+1.0, N2r-1.0, beta2 ) +
8726 ecoul( Z1UCD-1.0, N1r+1.0, beta1,
8727 Z2UCD+1.0, N2r-1.0, beta2, d );
8728 re2 =
frldm( Z1UCD, N1r, beta1) +
8729 frldm( Z2UCD, N2r, beta2 ) +
8730 ecoul( Z1UCD, N1r, beta1,
8731 Z2UCD, N2r, beta2, d );
8732 re3 =
frldm( Z1UCD+1.0, N1r-1.0, beta1 ) +
8733 frldm( Z2UCD-1.0, N2r+1.0, beta2 ) +
8734 ecoul( Z1UCD+1.0, N1r-1.0, beta1,
8735 Z2UCD-1.0, N2r+1.0, beta2, d );
8736 eps2 = ( re1 - 2.0*re2 + re3 ) / 2.0;
8737 eps1 = ( re3 - re1 ) / 2.0;
8738 DN1_POL = -eps1 / ( 2.0 * eps2 );
8740 Z1rr = Z1UCD + DN1_POL;
8745 DN1_POL = DN1_POL - 0.6 *
Uwash(E_scission_post,Ecrit,FREDSHELL,gamma);
8746 Z1rr = Z1UCD + DN1_POL;
8747 if ( Z1rr < 50. ) Z1rr = 50.0;
8749 DN1_POL = DN1_POL + 0.60 *
Uwash(E_scission_post,Ecrit,FREDSHELL,gamma);
8750 Z1rr = Z1UCD + DN1_POL;
8751 if ( Z1rr > 50.0 ) Z1rr = 50.0;
8761 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED" << std::endl;
8765 if( (
utilabs(Z1rr - Z1r) > 3.0*ZA1width) || (Z1r < 1.0) )
goto fiss2803;
8777 z2 =
dint( Z_scission ) -
z1;
8779 N2 =
dint( N_scission ) - N1;
8783 if( (z1 < 0) || (z2 < 0) || (a1 < 0) || (a2 < 0) ){
8784 std::cout <<
" -------------------------------" << std::endl;
8785 std::cout <<
" Z, A, N : " << Z <<
" " << A <<
" " << N << std::endl;
8786 std::cout << z1 <<
" " << z2 <<
" " << a1 <<
" " << a2 << std::endl;
8787 std::cout << E_scission_post <<
" " << A_levdens <<
" " << CZ << std::endl;
8789 std::cout <<
" -------------------------------" << std::endl;
8798 if( N1mean > N*0.50 ){
8802 if(beta2< beta2gs) beta2 = beta2gs;
8803 E1exc = E_scission_pre * a1 / A + E_defo;
8804 E_defo =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8805 E2exc = E_scission_pre * a2 / A + E_defo;
8809 if(beta1< beta1gs) beta1 = beta1gs;
8810 E_defo =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8811 E1exc = E_scission_pre * a1 / A + E_defo;
8813 E2exc = E_scission_pre * a2 / A + E_defo;
8820 if( N1mean > N*0.5 ){
8823 if(beta1< beta1gs) beta1 = beta1gs;
8824 E_defo =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8825 E1exc = E_scission_pre * a1 / A + E_defo;
8827 if(beta2< beta2gs) beta2 = beta2gs;
8828 E_defo =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8829 E2exc = E_scission_pre * a2 / A + E_defo;
8833 if(beta2< beta2gs) beta2 = beta2gs;
8834 E_defo =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8835 E2exc = E_scission_pre * a2 / A + E_defo;
8837 if(beta1< beta1gs) beta1 = beta1gs;
8838 E_defo =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8839 E1exc = E_scission_pre * a1 / A + E_defo;
8846 if(beta1< beta1gs) beta1 = beta1gs;
8848 if(beta2< beta2gs) beta2 = beta2gs;
8849 E_defo1 =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8850 E_defo2 =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8851 E1exc = E_scission_pre * a1 / A + E_defo1;
8852 E2exc = E_scission_pre * a2 / A + E_defo2;
8857 TKER = ( z1 * z2 * 1.440 ) /
8858 ( R0 * std::pow(a1,0.333330) * (1.0 + 2.0/3.0 * beta1 ) +
8859 R0 * std::pow(a2,0.333330) * (1.0 + 2.0/3.0 * beta2 ) + 2.0 );
8861 EkinR1 = TKER * a2 /
A;
8862 EkinR2 = TKER * a1 /
A;
8863 v1 = std::sqrt(EkinR1/a1) * 1.3887;
8864 v2 = std::sqrt(EkinR2/a2) * 1.3887;
8872 e1 =
gausshaz(0,E1exc,E1exc_sigma);
8873 if(e1<0.)
goto fis987;
8876 e2 =
gausshaz(0,E2exc,E2exc_sigma);
8877 if(e2<0.)
goto fis988;
8879 (*NbLam0_par) = NbLam0;
8917 G4double r_in = 0.0, r_rest = 0.0, r_help = 0.0;
8923 r_in = r_origin + 0.5;
8925 if (r_even_odd < 0.001) {
8926 i_out = (
G4int)(r_floor);
8929 r_rest = r_in - r_floor;
8930 r_middle = r_floor + 0.5;
8931 n_floor = (
G4int)(r_floor);
8932 if (n_floor%2 == 0) {
8934 r_help = r_middle + (r_rest - 0.5) * (1.0 - r_even_odd);
8938 r_help = r_middle + (r_rest - 0.5) * (1.0 + r_even_odd);
8940 i_out = (
G4int)(r_help);
8956 G4double xcom = 0.0, xvs = 0.0, xe = 0.0;
8960 alpha = ( std::sqrt(5.0/(4.0*pi)) ) * beta;
8962 xcom = 1.0 - 1.7826 * ((a - 2.0*
z)/a)*((a - 2.0*
z)/a);
8964 xvs = - xcom * ( 15.4941 * a -
8965 17.9439 * std::pow(a,2.0/3.0) * (1.0+0.4*alpha*
alpha) );
8967 xe = z*z * (0.7053/(std::pow(a,1.0/3.0)) * (1.0-0.2*alpha*alpha) - 1.1529/
a);
8994 dtot = r0 * ( std::pow((z1+n1),1.0/3.0) * (1.0+0.6666667*beta1)
8995 + std::pow((z2+n2),1.0/3.0) * (1.0+0.6666667*beta2) ) +
d;
8996 fecoul = z1 * z2 * 1.44 / dtot;
9008 R_wash = std::exp(-E * Freduction * gamma);
9010 R_wash = std::exp(- Ecrit * Freduction * gamma -(E-Ecrit) * gamma);
9068 G4double z = 0.0,
n = 0.0,
a = 0.0, av = 0.0, as = 0.0;
9069 G4double a0 = 0.0,
c1 = 0.0, c4 = 0.0, b1 = 0.0, b3 = 0.0;
9071 G4double r0 = 0.0, kf = 0.0, ks = 0.0;
9072 G4double kv = 0.0, rp = 0.0, ay = 0.0, aden = 0.0, x0 = 0.0, y0 = 0.0;
9073 G4double esq = 0.0, ael = 0.0, i = 0.0;
9124 c1 = 3.0/5.0*esq/r0;
9125 c4 = 5.0/4.0*std::pow((3.0/(2.0*pi)),(2.0/3.0)) *
c1;
9126 kf = std::pow((9.0*pi*z/(4.0*
a)),(1.0/3.0))/r0;
9128 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));
9132 x0 = r0 * std::pow(
a,(1.0/3.0)) / ay;
9133 y0 = r0 * std::pow(
a,(1.0/3.0)) / aden;
9135 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);
9137 b3 = 1.0 - 5.0/std::pow(y0,2) * (1.0 - 15.0/(8.0*y0) + 21.0/(8.0 * std::pow(y0,3))
9138 - 3.0/4.0 * (1.0 + 9.0/(2.0*y0) + 7.0/std::pow(y0,2)
9139 + 7.0/(2.0 * std::pow(y0,3))) * std::exp(-2.0*y0));
9143 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
9144 +
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))
9145 + ff*std::pow(z,2)/
a -ca*(
n-
z) - ael * std::pow(z,(2.39e0));
9151 return eflmacResult;
9156 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][6],
G4int *ILOOP){
9158 G4int INMIN,INMAX,NDIF=0,IMEM;
9159 G4int NEVA=0,PEVA=0;
9167 for(
G4int i=0;i<200;i++){
9168 BU_TAB_TEMP[i][0] = 0.0;
9169 BU_TAB_TEMP[i][1] = 0.0;
9170 BU_TAB_TEMP[i][2] = 0.0;
9171 BU_TAB_TEMP[i][3] = 0.0;
9172 BU_TAB_TEMP[i][4] = 0.0;
9179 if(AFP==0 && ZFP==0){
9183 if((AFP==1 && ZFP==0) ||
9184 (AFP==1 && ZFP==1) ||
9185 (AFP==2 && ZFP==1) ||
9186 (AFP==3 && ZFP==1) ||
9187 (AFP==3 && ZFP==2) ||
9188 (AFP==4 && ZFP==2) ||
9189 (AFP==6 && ZFP==2) ||
9198 if ((AFP-ZFP)==0 && ZFP>1){
9201 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9202 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9203 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9204 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9205 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9206 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9207 *ILOOP = *ILOOP + 1;
9242 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9243 BU_TAB_TEMP[
I+1+*ILOOP][0] = 1.0;
9244 BU_TAB_TEMP[
I+1+*ILOOP][1] = 1.0;
9245 BU_TAB_TEMP[
I+1+*ILOOP][2] = VP2X;
9246 BU_TAB_TEMP[
I+1+*ILOOP][3] = VP2Y;
9247 BU_TAB_TEMP[
I+1+*ILOOP][4] = VP2Z;
9252 *ILOOP = *ILOOP + IMEM;
9257 NEVA = NDIF - INMAX;
9266 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9268 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9269 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9270 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9271 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9272 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9273 *ILOOP = *ILOOP + 1;
9280 if ((AFP>=2) && (ZFP==0)){
9285 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9287 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9288 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9289 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9290 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9291 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9292 *ILOOP = *ILOOP + 1;
9304 std::cout <<
"WARNING - BU UNSTABLE: AF < ZF" << std::endl;
9309 if ((AFP>=4) && (ZFP==1)){
9315 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9317 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9318 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9319 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9320 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9321 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9322 *ILOOP = *ILOOP + 1;
9334 if ((AFP==4) && (ZFP==3)){
9342 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9344 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9345 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9346 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9347 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9348 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9349 *ILOOP = *ILOOP + 1;
9351 if ((AFP==5) && (ZFP==2)){
9359 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9360 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9361 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9362 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9363 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9364 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9365 *ILOOP = *ILOOP + 1;
9368 if ((AFP==5) && (ZFP==3)){
9376 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9377 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9378 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9379 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9380 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9381 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9382 *ILOOP = *ILOOP + 1;
9385 if ((AFP==6) && (ZFP==4)){
9394 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9395 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9396 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9397 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9398 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9399 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9400 *ILOOP = *ILOOP + 1;
9408 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9409 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9410 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9411 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9412 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9413 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9414 *ILOOP = *ILOOP + 1;
9416 if ((AFP==7)&&(ZFP==2)){
9424 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9425 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9426 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9427 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9428 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9429 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9430 *ILOOP = *ILOOP + 1;
9433 if ((AFP==7) && (ZFP==5)){
9435 for(
int I = 0;
I<= AFP-5;
I++){
9437 double(AFP-
I-1),
double(ZFP-
I-1),
9439 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9440 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9441 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9442 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9443 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9444 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9445 *ILOOP = *ILOOP + 1;
9456 if ((AFP==8) && (ZFP==4)){
9463 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9464 BU_TAB_TEMP[*ILOOP][0] = 2.0;
9465 BU_TAB_TEMP[*ILOOP][1] = 4.0;
9466 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9467 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9468 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9469 *ILOOP = *ILOOP + 1;
9471 if ((AFP==8) && (ZFP==6)){
9480 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9481 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9482 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9483 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9484 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9485 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9486 *ILOOP = *ILOOP + 1;
9493 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9494 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9495 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9496 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9497 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9498 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9499 *ILOOP = *ILOOP + 1;
9505 if((AFP==9) && (ZFP==2)){
9514 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9515 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9516 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9517 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9518 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9519 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9520 *ILOOP = *ILOOP + 1;
9526 if((AFP==9) && (ZFP==5)){
9534 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9535 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9536 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9537 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9538 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9539 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9540 *ILOOP = *ILOOP + 1;
9547 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9548 BU_TAB_TEMP[*ILOOP][0] = 2.0;
9549 BU_TAB_TEMP[*ILOOP][1] = 4.0;
9550 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9551 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9552 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9553 *ILOOP = *ILOOP + 1;
9559 if((AFP==10) && (ZFP==2)){
9568 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9569 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9570 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9571 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9572 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9573 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9574 *ILOOP = *ILOOP + 1;
9582 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9583 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9584 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9585 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9586 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9587 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9588 *ILOOP = *ILOOP + 1;
9593 if ((AFP==10) && (ZFP==3)){
9601 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9602 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9603 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9604 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9605 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9606 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9607 *ILOOP = *ILOOP + 1;
9612 if ((AFP==10) && (ZFP==7)){
9620 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9621 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9622 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9623 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9624 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9625 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9626 *ILOOP = *ILOOP + 1;
9632 if((AFP==11) && (ZFP==7)){
9640 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9641 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9642 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9643 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9644 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9645 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9646 *ILOOP = *ILOOP + 1;
9651 if ((AFP==12) && (ZFP==8)){
9660 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9661 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9662 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9663 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9664 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9665 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9666 *ILOOP = *ILOOP + 1;
9673 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9674 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9675 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9676 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9677 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9678 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9679 *ILOOP = *ILOOP + 1;
9684 if ((AFP==15) && (ZFP==9)){
9692 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9693 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9694 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9695 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9696 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9697 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9698 *ILOOP = *ILOOP + 1;
9704 if ((AFP==16) && (ZFP==9)){
9712 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9713 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9714 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9715 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9716 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9717 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9718 *ILOOP = *ILOOP + 1;
9724 if ((AFP==16) && (ZFP==10)){
9732 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9733 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9734 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9735 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9736 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9737 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9738 *ILOOP = *ILOOP + 1;
9745 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9746 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9747 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9748 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9749 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9750 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9751 *ILOOP = *ILOOP + 1;
9756 if((AFP==18) && (ZFP==11)){
9764 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9765 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9766 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9767 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9768 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9769 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9770 *ILOOP = *ILOOP + 1;
9775 if((AFP==19) && (ZFP==11)){
9783 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9784 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9785 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9786 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9787 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9788 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9789 *ILOOP = *ILOOP + 1;
9794 if (ZFP>=4 && (AFP-ZFP)==1){
9803 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9804 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9805 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9806 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9807 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9808 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9809 *ILOOP = *ILOOP + 1;
9818 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9819 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9820 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9821 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9822 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9823 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9824 *ILOOP = *ILOOP + 1;
9842 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){
9845 G4double PX1,PX2,PY1,PY2,PZ1,PZ2,PTOT;
9846 G4double RNDT,CTET1,STET1,RNDP,PHI1,ETOT_P1,ETOT_P2;
9848 G4double vxout=0.,vyout=0.,vzout=0.;
9849 G4int iain,izin,ianew,iznew,inin,innew;
9859 innew = ianew - iznew;
9862 mglms(ain,zin,3,&MASS);
9863 mglms(anew,znew,3,&MASS1);
9864 mglms(ain-anew,zin-znew,3,&MASS2);
9865 ekin_tot = MASS-MASS1-MASS2;
9869 if(izin>12)std::cout <<
"*** ZIN > 12 ***" << izin << std::endl;
9872 if( ekin_tot<0.00 ){
9883 EKIN_P1 = ekin_tot*(ain-anew)/ ain;
9884 ETOT_P1 = EKIN_P1 + anew * AMU;
9885 PTOT = anew*AMU*std::sqrt((EKIN_P1/(anew*AMU)+1.0)*(EKIN_P1/(anew*AMU)+1.0)-1.0);
9888 CTET1 = 2.0*RNDT-1.0;
9889 STET1 = std::sqrt(1.0-CTET1*CTET1);
9891 PHI1 = RNDP*2.0*3.141592654;
9892 PX1 = PTOT * STET1*std::cos(PHI1);
9893 PY1 = PTOT * STET1*std::sin(PHI1);
9895 *v1x = C * PX1 / ETOT_P1;
9896 *v1y = C * PY1 / ETOT_P1;
9897 *v1z = C * PZ1 / ETOT_P1;
9898 lorentz_boost(vxin,vyin,vzin,*v1x,*v1y,*v1z,&vxout,&vyout,&vzout);
9906 ETOT_P2 = (ekin_tot - EKIN_P1) + (ain-anew) * AMU;
9907 *v2x = C * PX2 / ETOT_P2;
9908 *v2y = C * PY2 / ETOT_P2;
9909 *v2z = C * PZ2 / ETOT_P2;
9910 lorentz_boost(vxin,vyin,vzin,*v2x,*v2y,*v2z,&vxout,&vyout,&vzout);
9928 G4double GAMMA,VR,
C,CC,DENO,VXNOM,VYNOM,VZNOM;
9939 VR = std::sqrt(VXR*VXR + VYR*VYR + VZR*VZR);
9946 GAMMA = 1.0/std::sqrt(1.0 - VR*VR/CC);
9947 DENO = 1.0 - VXR*VXIN/CC - VYR*VYIN/CC - VZR*VZIN/CC;
9950 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;
9952 *VXOUT = VXNOM / (GAMMA * DENO);
9955 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;
9957 *VYOUT = VYNOM / (GAMMA * DENO);
9960 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;
9962 *VZOUT = VZNOM / (GAMMA * DENO);
9974 G4double EFF1=0.,EFF2=0.,VFF1=0.,VFF2=0.,
9975 AF1=0.,ZF1=0.,AF2=0.,ZF2=0.,
9976 AFF1=0.,ZFF1=0.,AFF2=0.,ZFF2=0.,
9977 vz1_eva=0., vx1_eva=0.,vy1_eva=0.,
9978 vz2_eva=0., vx2_eva=0.,vy2_eva=0.,
9979 vx_eva_sc=0.,vy_eva_sc=0.,vz_eva_sc=0.,
9980 VXOUT=0.,VYOUT=0.,VZOUT=0.,
9981 VX2OUT=0.,VY2OUT=0.,VZ2OUT=0.;
9982 G4int IEV_TAB_FIS=0,IEV_TAB_TEMP=0;
9983 G4double EV_TEMP1[200][6], EV_TEMP2[200][6],mtota=0.;
9984 G4int inttype = 0,inum=0;
9987 G4int NbLam0= (*NbLam0_par);
9989 for(
G4int I1=0;I1<200;I1++)
9990 for(
G4int I2=0;I2<6;I2++){
9991 EV_TEMP[I1][I2] = 0.0;
9992 EV_TEMP1[I1][I2] = 0.0;
9993 EV_TEMP2[I1][I2] = 0.0;
9996 G4double et = EE - JPRF * JPRF * 197. * 197./(2.*0.4*931.*std::pow(AF,5.0/3.0)*1.16*1.16);
9998 fissionDistri(AF,ZF,et,AF1,ZF1,EFF1,VFF1,AF2,ZF2,EFF2,VFF2,
9999 vx_eva_sc,vy_eva_sc,vz_eva_sc,&NbLam0);
10004 G4double pbH = (AF1 - ZF1) / (AF1 - ZF1 + AF2 - ZF2);
10005 for(
G4int i=0;i<NbLam0;i++){
10025 G4double VPERP1 = std::sqrt(VFF1*VFF1 - VZ1_FISSION*VZ1_FISSION);
10027 G4double VX1_FISSION = VPERP1 * std::sin(ALPHA1);
10028 G4double VY1_FISSION = VPERP1 * std::cos(ALPHA1);
10029 G4double VX2_FISSION = - VX1_FISSION / VFF1 * VFF2;
10030 G4double VY2_FISSION = - VY1_FISSION / VFF1 * VFF2;
10031 G4double VZ2_FISSION = - VZ1_FISSION / VFF1 * VFF2;
10034 if( (ZF1<=0.0) || (AF1<=0.0) || (AF1<ZF1) ){
10035 std::cout <<
"F1 unphysical: "<<ZF<<
" "<<AF<<
" "<<EE<<
" "<<ZF1<<
" "<<AF1 << std::endl;
10041 G4int FF11=0, FIMF11=0;
10042 G4double ZIMFF1=0., AIMFF1=0.,TKEIMF1=0.,JPRFOUT=0.;
10044 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,&NbLam1);
10046 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
10047 EV_TEMP[IJ+IEV_TAB_FIS][0] = EV_TEMP1[IJ][0];
10048 EV_TEMP[IJ+IEV_TAB_FIS][1] = EV_TEMP1[IJ][1];
10055 EV_TEMP1[IJ][2],EV_TEMP1[IJ][3],EV_TEMP1[IJ][4],
10056 &VXOUT,&VYOUT,&VZOUT);
10059 &VX2OUT,&VY2OUT,&VZ2OUT);
10060 EV_TEMP[IJ+IEV_TAB_FIS][2] = VX2OUT;
10061 EV_TEMP[IJ+IEV_TAB_FIS][3] = VY2OUT;
10062 EV_TEMP[IJ+IEV_TAB_FIS][4] = VZ2OUT;
10065 IEV_TAB_FIS = IEV_TAB_FIS + IEV_TAB_TEMP;
10070 if( (ZF2<=0.0) || (AF2<=0.0) || (AF2<ZF2) ){
10071 std::cout <<
"F2 unphysical: "<<ZF<<
" "<<AF<<
" "<<EE<<
" "<<ZF2<<
" "<<AF2 << std::endl;
10077 G4int FF22=0, FIMF22=0;
10078 G4double ZIMFF2=0., AIMFF2=0.,TKEIMF2=0.,JPRFOUT=0.;
10080 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,&NbLam2);
10082 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
10083 EV_TEMP[IJ+IEV_TAB_FIS][0] = EV_TEMP2[IJ][0];
10084 EV_TEMP[IJ+IEV_TAB_FIS][1] = EV_TEMP2[IJ][1];
10091 EV_TEMP2[IJ][2],EV_TEMP2[IJ][3],EV_TEMP2[IJ][4],
10092 &VXOUT,&VYOUT,&VZOUT);
10095 &VX2OUT,&VY2OUT,&VZ2OUT);
10096 EV_TEMP[IJ+IEV_TAB_FIS][2] = VX2OUT;
10097 EV_TEMP[IJ+IEV_TAB_FIS][3] = VY2OUT;
10098 EV_TEMP[IJ+IEV_TAB_FIS][4] = VZ2OUT;
10101 IEV_TAB_FIS = IEV_TAB_FIS + IEV_TAB_TEMP;
10114 VX1_FISSION,VY1_FISSION,VZ1_FISSION,
10115 &VXOUT,&VYOUT,&VZOUT);
10116 VX1_FISSION = VXOUT;
10117 VY1_FISSION = VYOUT;
10118 VZ1_FISSION = VZOUT;
10120 VX2_FISSION,VY2_FISSION,VZ2_FISSION,
10121 &VXOUT,&VYOUT,&VZOUT);
10122 VX2_FISSION = VXOUT;
10123 VY2_FISSION = VYOUT;
10124 VZ2_FISSION = VZOUT;
10129 (*VX1_FISSION_par) = VX1_FISSION;
10130 (*VY1_FISSION_par) = VY1_FISSION;
10131 (*VZ1_FISSION_par) = VZ1_FISSION;
10132 (*VX_EVA_SC_par)=vx_eva_sc;
10133 (*VY_EVA_SC_par)=vy_eva_sc;
10134 (*VZ_EVA_SC_par)=vz_eva_sc;
10138 (*VX2_FISSION_par) = VX2_FISSION;
10139 (*VY2_FISSION_par) = VY2_FISSION;
10140 (*VZ2_FISSION_par) = VZ2_FISSION;
10141 (*IEV_TAB_FIS_par) = IEV_TAB_FIS;
10142 (*NbLam0_par) = NbLam1 + NbLam2;
10143 if(NbLam0>(NbLam1 + NbLam2))
varntp->
kfis = 25;
10150 G4double V_over_V0,R0,RALL,RHAZ,
R,TKE,Ekin,V,VPERP,ALPHA1;
10162 RALL = R0 * std::pow(V_over_V0,1.0/3.0) * std::pow(AAL,1.0/3.0);
10164 R = std::pow(RHAZ,1.0/3.0) * RALL;
10165 TKE = 1.44 * Z * ZALL * R*R * (1.0 - A/AAL)*(1.0 - A/AAL) / std::pow(RALL,3.0);
10167 Ekin = TKE * (AAL -
A) / AAL;
10169 V = std::sqrt(Ekin/A) * 1.3887;
10171 VPERP = std::sqrt(V*V - (*VZ)*(*VZ));
10173 *VX = VPERP * std::sin(ALPHA1);
10174 *VY = VPERP * std::cos(ALPHA1);
10200 ix =
G4int(
y * 100 + 43543000);
10201 if(
mod(ix,2) == 0) {
10215 l_plus = lambda + 1.;
10219 y=std::pow(
G4AblaRandom::flat()*(std::pow(rxmax,l_plus)-std::pow(rxmin,l_plus))+ std::pow(rxmin,l_plus),1.0/l_plus);
10253 G4double Efermi = 5.0 * SIGMA_0 * SIGMA_0 / (2.0 * 931.4940);
10259 SIGMA_0 = SIGMA_0 * std::pow(V0_over_VBU,1.0/3.0);
10264 GOLDHA_BU = SIGMA_0 * std::sqrt((APRF*(AABRA-APRF))/(AABRA-1.0));
10265 GOLDHA = GOLDHA_BU*std::sqrt(1.0 +
10273 (AABRA - APRF) / AABRA);
10274 GOLDHA = GOLDHA_BU;
10279 GOLDHA = SIGMA_0 * std::sqrt((APRF*(AABRA-APRF))/(AABRA-1.0));
10287 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING PX IN Rn07.FOR. A VALUE WILL BE FORCED." << std::endl;
10288 *PX = (AABRA-1.0)*931.4940;
10290 if(
std::abs(*PX)>= AABRA*931.494){
10299 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING PY IN Rn07.FOR. A VALUE WILL BE FORCED." << std::endl;
10300 *PY = (AABRA-1.0)*931.4940;
10302 if(
std::abs(*PY)>= AABRA*931.494){
10311 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING PZ IN Rn07.FOR. A VALUE WILL BE FORCED." << std::endl;
10312 *PZ = (AABRA-1.0)*931.4940;
10314 if(
std::abs(*PZ)>= AABRA*931.494){
10331 v1 = 2.0*
haz(k) - 1.0;
10332 v2 = 2.0*
haz(k) - 1.0;
10333 r = std::pow(v1,2) + std::pow(v2,2);
10336 fac = std::sqrt(-2.*std::log(r)/r);
10338 fgausshaz = v2*fac*sig+xmoy;
10342 fgausshaz=gset*sig+xmoy;