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;
115 G4int NbLam0 = std::abs(nucleusS);
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++){
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;
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];
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);
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;
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 /
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;
884 if(std::abs(PX_BU_SUM)>10. || std::abs(PY_BU_SUM)>10. ||
885 std::abs(PZ_BU_SUM)>10.){
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 /
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;
1105 G4double VX1_IMF = VPERP1 * std::sin(ALPHA1);
1106 G4double VY1_IMF = VPERP1 * std::cos(ALPHA1);
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];
1257 &VP1X,&VP1Y,&VP1Z,BU_TAB_TEMP1,&ILOOP);
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;
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);
1320 if(IMULTBU>200)std::cout <<
"IMULTBU>200 " << IMULTBU << std::endl;
1351 if(zprfp<=2 && zprfp<aprfp){
1369 if(zprfp<=2 && zprfp==aprfp){
1371 VX_PREF, VY_PREF, VZ_PREF,
1372 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
1378 for(
G4int I = 0;I<ILOOP;I++){
1379 for(
G4int IJ = 0; IJ<6; IJ++)
1380 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1382 IEV_TAB = IEV_TAB + ILOOP;
1400 VX_PREF, VY_PREF, VZ_PREF,
1401 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
1407 for(
G4int I = 0;I<ILOOP;I++){
1408 for(
G4int IJ = 0; IJ<6; IJ++)
1409 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1411 IEV_TAB = IEV_TAB + ILOOP;
1426 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);
1428 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1429 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1430 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1431 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1438 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1439 &VXOUT,&VYOUT,&VZOUT);
1440 EV_TAB[IJ+IEV_TAB][2] = VXOUT;
1441 EV_TAB[IJ+IEV_TAB][3] = VYOUT;
1442 EV_TAB[IJ+IEV_TAB][4] = VZOUT;
1444 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1450 vx_eva,vy_eva,vz_eva,
1451 &VXOUT,&VYOUT,&VZOUT);
1456 if(ff == 0 && fimf == 0){
1468 VFP1_CM[0] = V_CM[0];
1469 VFP1_CM[1] = V_CM[1];
1470 VFP1_CM[2] = V_CM[2];
1471 for(
G4int j=0;j<3;j++){
1477 if(ff == 1 && fimf == 0) ftype = 1;
1478 if(ff == 0 && fimf == 1) ftype = 2;
1491 G4int IEV_TAB_FIS = 0,imode=0;
1493 G4double vx1_fission=0.,vy1_fission=0.,vz1_fission=0.;
1494 G4double vx2_fission=0.,vy2_fission=0.,vz2_fission=0.;
1495 G4double vx_eva_sc=0.,vy_eva_sc=0.,vz_eva_sc=0.;
1498 &vx1_fission,&vy1_fission,&vz1_fission,
1499 &vx2_fission,&vy2_fission,&vz2_fission,
1500 &ZFP1,&AFP1,&SFP1,&ZFP2,&AFP2,&SFP2,&imode,
1501 &vx_eva_sc,&vy_eva_sc,&vz_eva_sc,EV_TEMP,&IEV_TAB_FIS,&NbLam0);
1503 for(
G4int IJ = 0; IJ< IEV_TAB_FIS;IJ++){
1504 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1505 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1506 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1513 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1514 &VXOUT,&VYOUT,&VZOUT);
1515 EV_TAB[IJ+IEV_TAB][2] = VXOUT;
1516 EV_TAB[IJ+IEV_TAB][3] = VYOUT;
1517 EV_TAB[IJ+IEV_TAB][4] = VZOUT;
1519 IEV_TAB = IEV_TAB + IEV_TAB_FIS;
1535 V_CM[0],V_CM[1],V_CM[2],
1536 &VXOUT,&VYOUT,&VZOUT);
1539 &VX2OUT,&VY2OUT,&VZ2OUT);
1540 VFP1_CM[0] = VX2OUT;
1541 VFP1_CM[1] = VY2OUT;
1542 VFP1_CM[2] = VZ2OUT;
1549 V_CM[0],V_CM[1],V_CM[2],
1550 &VXOUT,&VYOUT,&VZOUT);
1553 &VX2OUT,&VY2OUT,&VZ2OUT);
1554 VFP2_CM[0] = VX2OUT;
1555 VFP2_CM[1] = VY2OUT;
1556 VFP2_CM[2] = VZ2OUT;
1560 }
else if(ftype == 2){
1569 G4double pbH = (af-zf) / (af-zf+aimf-zimf);
1571 for(
G4int i=0;i<NbLam0;i++){
1580 G4double EkinR1 = tkeimf * aimf / (af+aimf);
1581 G4double EkinR2 = tkeimf * af / (af+aimf);
1583 G4double V2 = std::sqrt(EkinR2/aimf) * 1.3887;
1587 G4double VX1_IMF = VPERP1 * std::sin(ALPHA1);
1588 G4double VY1_IMF = VPERP1 * std::cos(ALPHA1);
1593 G4double EEIMFP = ee * af /(af + aimf);
1594 G4double EEIMF = ee * aimf /(af + aimf);
1597 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.));
1599 G4double JPRFHEAVY = jprf0 * 0.4 * 931.49 * 1.16*1.16 * std::pow(af,5.0/3.0) / IINERTTOT;
1600 G4double JPRFLIGHT = jprf0 * 0.4 * 931.49 * 1.16*1.16 * std::pow(aimf,5.0/3.0) / IINERTTOT;
1601 if(af<2.0) std::cout <<
"RN117-4,AF,ZF,EE,JPRFheavy" << std::endl;
1603 G4double vx1ev_imf=0., vy1ev_imf=0., vz1ev_imf=0., zdummy=0., adummy=0., tkedummy=0.,jprf1=0.;
1605 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);
1607 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1608 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1609 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1610 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1617 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1618 &VXOUT,&VYOUT,&VZOUT);
1621 &VX2OUT,&VY2OUT,&VZ2OUT);
1622 EV_TAB[IJ+IEV_TAB][2] = VX2OUT;
1623 EV_TAB[IJ+IEV_TAB][3] = VY2OUT;
1624 EV_TAB[IJ+IEV_TAB][4] = VZ2OUT;
1626 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1635 G4double zffimf, affimf,zdummy1=0., adummy1=0., tkedummy1=0.,jprf2,vx2ev_imf,vy2ev_imf,
1638 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);
1640 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1641 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1642 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1643 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1650 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1651 &VXOUT,&VYOUT,&VZOUT);
1654 &VX2OUT,&VY2OUT,&VZ2OUT);
1655 EV_TAB[IJ+IEV_TAB][2] = VX2OUT;
1656 EV_TAB[IJ+IEV_TAB][3] = VY2OUT;
1657 EV_TAB[IJ+IEV_TAB][4] = VZ2OUT;
1659 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1673 V_CM[0],V_CM[1],V_CM[2],
1674 &VXOUT,&VYOUT,&VZOUT);
1677 &VX2OUT,&VY2OUT,&VZ2OUT);
1678 VIMF_CM[0] = VX2OUT;
1679 VIMF_CM[1] = VY2OUT;
1680 VIMF_CM[2] = VZ2OUT;
1686 V_CM[0],V_CM[1],V_CM[2],
1687 &VXOUT,&VYOUT,&VZOUT);
1690 &VX2OUT,&VY2OUT,&VZ2OUT);
1691 VFP1_CM[0] = VX2OUT;
1692 VFP1_CM[1] = VY2OUT;
1693 VFP1_CM[2] = VZ2OUT;
1695 if(FF11==0 && FIMF11==0){
1707 for(
G4int I=0;I<3;I++)
1711 }
else if(FF11==1 && FIMF11==0){
1725 G4int IEV_TAB_FIS = 0,imode=0;
1727 G4double vx1_fission=0.,vy1_fission=0.,vz1_fission=0.;
1728 G4double vx2_fission=0.,vy2_fission=0.,vz2_fission=0.;
1729 G4double vx_eva_sc=0.,vy_eva_sc=0.,vz_eva_sc=0.;
1732 &vx1_fission,&vy1_fission,&vz1_fission,
1733 &vx2_fission,&vy2_fission,&vz2_fission,
1734 &ZFP1,&AFP1,&SFP1,&ZFP2,&AFP2,&SFP2,&imode,
1735 &vx_eva_sc,&vy_eva_sc,&vz_eva_sc,EV_TEMP,&IEV_TAB_FIS,&NbLamH);
1737 for(
int IJ = 0; IJ< IEV_TAB_FIS;IJ++){
1738 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1739 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1740 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1747 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1748 &VXOUT,&VYOUT,&VZOUT);
1749 EV_TAB[IJ+IEV_TAB][2] = VXOUT;
1750 EV_TAB[IJ+IEV_TAB][3] = VYOUT;
1751 EV_TAB[IJ+IEV_TAB][4] = VZOUT;
1753 IEV_TAB = IEV_TAB + IEV_TAB_FIS;
1766 V_CM[0],V_CM[1],V_CM[2],
1767 &VXOUT,&VYOUT,&VZOUT);
1770 &VX2OUT,&VY2OUT,&VZ2OUT);
1772 VX2OUT,VY2OUT,VZ2OUT,
1773 &VXOUT,&VYOUT,&VZOUT);
1776 &VX2OUT,&VY2OUT,&VZ2OUT);
1777 VFP1_CM[0] = VX2OUT;
1778 VFP1_CM[1] = VY2OUT;
1779 VFP1_CM[2] = VZ2OUT;
1789 V_CM[0],V_CM[1],V_CM[2],
1790 &VXOUT,&VYOUT,&VZOUT);
1793 &VX2OUT,&VY2OUT,&VZ2OUT);
1795 VX2OUT,VY2OUT,VZ2OUT,
1796 &VXOUT,&VYOUT,&VZOUT);
1799 &VX2OUT,&VY2OUT,&VZ2OUT);
1800 VFP2_CM[0] = VX2OUT;
1801 VFP2_CM[1] = VY2OUT;
1802 VFP2_CM[2] = VZ2OUT;
1804 }
else if(FF11==0 && FIMF11==1){
1821 G4double pbH1 = (af-zf) / (af-zf+aimf-zimf);
1822 for(
G4int i=0;i<NbLamH;i++){
1831 EkinR1 = tkeimf * aimf / (af+aimf);
1832 EkinR2 = tkeimf * af / (af+aimf);
1833 V1 = std::sqrt(EkinR1/af) * 1.3887;
1834 V2 = std::sqrt(EkinR2/aimf) * 1.3887;
1836 VPERP1 = std::sqrt(
V1*
V1 - VZ1_IMFS*VZ1_IMFS);
1838 G4double VX1_IMFS = VPERP1 * std::sin(ALPHA1);
1839 G4double VY1_IMFS = VPERP1 * std::cos(ALPHA1);
1844 EEIMFP = ee * af /(af + aimf);
1845 EEIMF = ee * aimf /(af + aimf);
1848 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.));
1850 JPRFHEAVY = jprf1 * 0.4 * 931.49 * 1.16*1.16 * std::pow(af,5.0/3.0) / IINERTTOT;
1851 JPRFLIGHT = jprf1 * 0.4 * 931.49 * 1.16*1.16 * std::pow(aimf,5.0/3.0) / IINERTTOT;
1853 G4double zffs=0.,affs=0.,vx1ev_imfs=0.,vy1ev_imfs=0.,vz1ev_imfs=0.,jprf3=0.;
1855 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);
1857 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1858 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1859 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1860 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1867 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1868 &VXOUT,&VYOUT,&VZOUT);
1871 &VX2OUT,&VY2OUT,&VZ2OUT);
1872 EV_TAB[IJ+IEV_TAB][2] = VX2OUT;
1873 EV_TAB[IJ+IEV_TAB][3] = VY2OUT;
1874 EV_TAB[IJ+IEV_TAB][4] = VZ2OUT;
1876 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1885 G4double zffimfs=0.,affimfs=0.,vx2ev_imfs=0.,vy2ev_imfs=0.,vz2ev_imfs=0.,jprf4=0.;
1887 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);
1889 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
1890 EV_TAB[IJ+IEV_TAB][0] = EV_TEMP[IJ][0];
1891 EV_TAB[IJ+IEV_TAB][1] = EV_TEMP[IJ][1];
1892 EV_TAB[IJ+IEV_TAB][5] = EV_TEMP[IJ][5];
1899 EV_TEMP[IJ][2],EV_TEMP[IJ][3],EV_TEMP[IJ][4],
1900 &VXOUT,&VYOUT,&VZOUT);
1903 &VX2OUT,&VY2OUT,&VZ2OUT);
1904 EV_TAB[IJ+IEV_TAB][2] = VX2OUT;
1905 EV_TAB[IJ+IEV_TAB][3] = VY2OUT;
1906 EV_TAB[IJ+IEV_TAB][4] = VZ2OUT;
1908 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1923 V_CM[0],V_CM[1],V_CM[2],
1924 &VXOUT,&VYOUT,&VZOUT);
1927 &VX2OUT,&VY2OUT,&VZ2OUT);
1929 VX2OUT,VY2OUT,VZ2OUT,
1930 &VXOUT,&VYOUT,&VZOUT);
1933 &VX2OUT,&VY2OUT,&VZ2OUT);
1934 VFP1_CM[0] = VX2OUT;
1935 VFP1_CM[1] = VY2OUT;
1936 VFP1_CM[2] = VZ2OUT;
1944 V_CM[0],V_CM[1],V_CM[2],
1945 &VXOUT,&VYOUT,&VZOUT);
1948 &VX2OUT,&VY2OUT,&VZ2OUT);
1950 VX2OUT,VY2OUT,VZ2OUT,
1951 &VXOUT,&VYOUT,&VZOUT);
1954 &VX2OUT,&VY2OUT,&VZ2OUT);
1955 VFP2_CM[0] = VX2OUT;
1956 VFP2_CM[1] = VY2OUT;
1957 VFP2_CM[2] = VZ2OUT;
1962 if(ftype!=1 && ftype!=21){
1968 VFP1_CM[0],VFP1_CM[1],VFP1_CM[2],
1969 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
1977 for(
G4int I = 0;I<ILOOP;I++){
1978 for(
G4int IJ = 0; IJ<5; IJ++)
1979 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1981 IEV_TAB = IEV_TAB + ILOOP;
1988 VIMF_CM[0],VIMF_CM[1],VIMF_CM[2],
1989 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
1997 for(
G4int I = 0;I<ILOOP;I++){
1998 for(
G4int IJ = 0; IJ<5; IJ++)
1999 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2001 IEV_TAB = IEV_TAB + ILOOP;
2008 VFP2_CM[0],VFP2_CM[1],VFP2_CM[2],
2009 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
2017 for(
G4int I = 0;I<ILOOP;I++){
2018 for(
G4int IJ = 0; IJ<5; IJ++)
2019 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2021 IEV_TAB = IEV_TAB + ILOOP;
2029 if(ftype==1 || ftype==21){
2034 VFP1_CM[0],VFP1_CM[1],VFP1_CM[2],
2035 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
2043 for(
G4int I = 0;I<ILOOP;I++){
2044 for(
G4int IJ = 0; IJ<5; IJ++)
2045 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2047 IEV_TAB = IEV_TAB + ILOOP;
2053 VFP2_CM[0],VFP2_CM[1],VFP2_CM[2],
2054 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
2062 for(
G4int I = 0;I<ILOOP;I++){
2063 for(
G4int IJ = 0; IJ<5; IJ++)
2064 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2066 IEV_TAB = IEV_TAB + ILOOP;
2073 VIMF_CM[0],VIMF_CM[1],VIMF_CM[2],
2074 &VP1X,&VP1Y,&VP1Z,EV_TAB_TEMP,&ILOOP);
2082 for(
G4int I = 0;I<ILOOP;I++){
2083 for(
G4int IJ = 0; IJ<5; IJ++)
2084 EV_TAB[I+IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2086 IEV_TAB = IEV_TAB + ILOOP;
2092 if((ftype == 1 || ftype == 21) && (AFP2<=0 || AFP1<=0 || ZFP2<=0 || ZFP1<=0)){
2093 std::cout <<
"ZFP1:" << ZFP1 << std::endl;
2094 std::cout <<
"AFP1:" << AFP1 << std::endl;
2095 std::cout <<
"ZFP2:" << ZFP2 << std::endl;
2096 std::cout <<
"AFP2:" << AFP2 << std::endl;
2100 EV_TAB[IEV_TAB][0] = ZFP1;
2101 EV_TAB[IEV_TAB][1] = AFP1;
2102 EV_TAB[IEV_TAB][5] = SFP1;
2103 EV_TAB[IEV_TAB][2] = VFP1_CM[0];
2104 EV_TAB[IEV_TAB][3] = VFP1_CM[1];
2105 EV_TAB[IEV_TAB][4] = VFP1_CM[2];
2106 IEV_TAB = IEV_TAB + 1;
2109 EV_TAB[IEV_TAB][0] = ZFP2;
2110 EV_TAB[IEV_TAB][1] = AFP2;
2111 EV_TAB[IEV_TAB][5] = SFP2;
2112 EV_TAB[IEV_TAB][2] = VFP2_CM[0];
2113 EV_TAB[IEV_TAB][3] = VFP2_CM[1];
2114 EV_TAB[IEV_TAB][4] = VFP2_CM[2];
2115 IEV_TAB = IEV_TAB + 1;
2119 EV_TAB[IEV_TAB][0] = ZFPIMF;
2120 EV_TAB[IEV_TAB][1] = AFPIMF;
2121 EV_TAB[IEV_TAB][5] = SFPIMF;
2122 EV_TAB[IEV_TAB][2] = VIMF_CM[0];
2123 EV_TAB[IEV_TAB][3] = VIMF_CM[1];
2124 EV_TAB[IEV_TAB][4] = VIMF_CM[2];
2125 IEV_TAB = IEV_TAB + 1;
2189#ifdef ABLAXX_IN_GEANT4_MODE
2194 if(dataInterface->
readData() ==
true) {
2203 for(
G4int z = 0; z < 99; z++) {
2214 for(
G4int z = 0; z < 137; z++){
2221 for(
G4int z = 0; z < 500; z++) {
2222 for(
G4int a = 0; a < 500; a++) {
2232 for(
G4int i=1;i<13;i++){
2233 for(
G4int j=1;j<154;j++){
2243 mfrldm[j][i] = MP*i+MN*j+
eflmac(i+j,i,1,0);
2248 for(
G4int i=1;i<13;i++){
2249 for(
G4int j=1;j<154;j++){
2269 e0 = 0.285+11.17*std::pow(j+i,-0.464) -0.390-0.00058*(j+i);
2275 e0 = 22.34*std::pow(j+i,-0.464)-0.235;
2282 if((j==i)&&
mod(j,2)==1&&
mod(i,2)==1){
2283 e0 = e0 - 30.0*(1.0/
G4double(j+i));
2297 delete dataInterface;
2392 G4double xv = 0.0, xs = 0.0, xc = 0.0, xa = 0.0;
2394 if ((a <= 0.01) || (z < 0.01)) {
2399 xs = 17.23*std::pow(a,(2.0/3.0));
2402 xc = 0.7*z*(z-1.0)*std::pow((a-1.0),(-1.e0/3.e0));
2409 xa = 23.6*(std::pow((a-2.0*z),2)/a);
2410 (*el) = xv+xs+xc+xa;
2438 if ( (a1 <= 0) || (z1 <= 0) || ((a1-z1) <= 0) ) {
2447 (*el) =
eflmac(a1,z1,0,refopt4);
2457 (*el) = (*el) + (12.552-0.1436*z1);
2459 if(n1>145&&n1<=152){
2460 (*el) = (*el) + ((152.4-1.77*z1)+(-0.972+0.0113*z1)*n1);
2480 G4double x = 0.0, v = 0.0, dx = 0.0;
2482 const G4int alpha2Size = 37;
2484 G4double alpha2[alpha2Size] = {0.0, 2.5464e0, 2.4944e0, 2.4410e0, 2.3915e0, 2.3482e0,
2485 2.3014e0, 2.2479e0, 2.1982e0, 2.1432e0, 2.0807e0, 2.0142e0, 1.9419e0,
2486 1.8714e0, 1.8010e0, 1.7272e0, 1.6473e0, 1.5601e0, 1.4526e0, 1.3164e0,
2487 1.1391e0, 0.9662e0, 0.8295e0, 0.7231e0, 0.6360e0, 0.5615e0, 0.4953e0,
2488 0.4354e0, 0.3799e0, 0.3274e0, 0.2779e0, 0.2298e0, 0.1827e0, 0.1373e0,
2489 0.0901e0, 0.0430e0, 0.0000e0};
2494 v = (x - 0.3)/dx + 1.0;
2505 return(
alpha2[index] + (
alpha2[index+1] -
alpha2[index]) / dx * ( x - (0.3e0 + dx*(index-1))));
2520 G4double aa = 0.0, zz = 0.0, i = 0.0,z2a,C_S,R,W,G,G1,G2,A_CC;
2526 z2a= zz*zz/aa-ny*(1115.-939.+sn-slam)/(0.7053*std::pow(a,2./3.));
2530 fissilityResult = std::pow(zz,2) / aa /50.8830e0 / (1.0e0 - 1.7826e0 * std::pow(i,2));
2535 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));
2540 fissilityResult = std::pow(zz,2) / aa /(48.e0*(1.e0 - 17.22e0*std::pow(i,4)));
2545 C_S = 21.13 * (1.0 - 2.3*i*i);
2546 R = 1.16 * std::pow(aa,1.0/3.0);
2548 G1 = 1.0 - 15.0/8.0*W+21.0/8.0*W*W*W;
2549 G2 = 1.0 + 9.0/2.0*W + 7.0*W*W + 7.0/2.0*W*W*W;
2550 G = 1.0 - 5.0*W*W*(G1 - 3.0/4.0*G2*std::exp(-2.0/W));
2551 A_CC = 3.0/5.0 * 1.44 * G / 1.16;
2552 fissilityResult = z2a * A_CC/(2.0*C_S);
2555 if (fissilityResult > 1.0) {
2556 fissilityResult = 1.0;
2559 if (fissilityResult < 0.0) {
2560 fissilityResult = 0.0;
2563 return fissilityResult;
2566void 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,
2567G4int *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)
2577 G4int ff = (*ff_par);
2578 G4int fimf = (*fimf_par);
2580 G4int inttype = (*inttype_par);
2581 G4int inum = (*inum_par);
2582 G4int NbLam0 = (*NbLam0_par);
2674 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;
2675 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;
2676 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;
2677 G4double zimf= 0.0,aimf= 0.0,bimf= 0.0,sbimf= 0.0,timf= 0.0;
2678 G4int itest = 0, sortie=0;
2680 G4double ctet1 = 0.0, stet1 = 0.0, phi1 = 0.0;
2684 G4int fgamma = 0, gammadecay = 0, flamb0decay=0;
2686 G4double jprfn=0.0, jprfp=0.0, jprfd=0.0, jprft=0.0, jprfhe=0.0, jprfa=0.0, jprflamb0 = 0.0;
2692 const G4double mu2 = 931.494*931.494;
2697 G4int IEV_TAB_TEMP=0;
2699 for(
G4int I1=0;I1<200;I1++)
2700 for(
G4int I2=0;I2<6;I2++)
2701 EV_TEMP[I1][I2] = 0.0;
2711 if(ee<0.|| zf<3.)
goto evapora100;
2712 direct(zf,af,ee,jprf,&probp,&probd,&probt,&probn,&probhe,&proba,&probg,&probimf,&probf,&problamb0,&ptotl,
2713 &sn,&sbp,&sbd,&sbt,&sbhe,&sba,&slamb0,
2714 &ecn,&ecp,&ecd,&ect,&eche,&eca,&ecg,&eclamb0,
2715 &
bp,&bd,&bt,&bhe,&ba,&
sp,&sd,&st,&she,&sa,&ef,&ts1,inttype,inum,itest,&sortie,&tcn,
2716 &jprfn, &jprfp, &jprfd, &jprft, &jprfhe, &jprfa, &jprflamb0, &tsum, NbLam0);
2720 if(ptotl==0.0)
goto evapora100;
2724 if(e>1e30)std::cout <<
"ERROR AT THE EXIT OF EVAPORA,E>1.D30,AF="<< af <<
" ZF=" << zf << std::endl;
2731 pc = std::sqrt(std::pow((1.0 + (ecn)/9.3956e2),2.) - 1.0) * 9.3956e2;
2738 else if(probp!=0.0){
2742 pc = std::sqrt(std::pow((1.0 + ecp/9.3827e2),2.) - 1.0) * 9.3827e2;
2749 else if(probd!=0.0){
2753 pc = std::sqrt(std::pow((1.0 + ecd/1.875358e3),2) - 1.0) * 1.875358e3;
2760 else if(probt!=0.0){
2764 pc = std::sqrt(std::pow((1.0 + ect/2.80828e3),2) - 1.0) * 2.80828e3;
2771 else if(probhe!=0.0){
2774 epsiln = she + eche;
2775 pc = std::sqrt(std::pow((1.0 + eche/2.80826e3),2) - 1.0) * 2.80826e3;
2782 else{
if(proba!=0.0){
2786 pc = std::sqrt(std::pow((1.0 + eca/3.72834e3),2) - 1.0) * 3.72834e3;
2808 pc = std::sqrt(std::pow((1.0 + eca/3.72834e3),2) - 1.0) * 3.72834e3;
2817 else if (x < proba+probhe) {
2821 epsiln = she + eche;
2822 pc = std::sqrt(std::pow((1.0 + eche/2.80826e3),2) - 1.0) * 2.80826e3;
2831 else if (x < proba+probhe+probt) {
2836 pc = std::sqrt(std::pow((1.0 + ect/2.80828e3),2) - 1.0) * 2.80828e3;
2845 else if (x < proba+probhe+probt+probd) {
2850 pc = std::sqrt(std::pow((1.0 + ecd/1.875358e3),2) - 1.0) * 1.875358e3;
2859 else if (x < proba+probhe+probt+probd+probp) {
2864 pc = std::sqrt(std::pow((1.0 + ecp/9.3827e2),2) - 1.0) * 9.3827e2;
2873 else if (x < proba+probhe+probt+probd+probp+probn) {
2878 pc = std::sqrt(std::pow((1.0 + (ecn)/9.3956e2),2.) - 1.0) * 9.3956e2;
2887 else if (x < proba+probhe+probt+probd+probp+probn+problamb0) {
2891 epsiln = slamb0 + eclamb0;
2892 pc = std::sqrt(std::pow((1.0 + (eclamb0)/11.1568e2),2.) - 1.0) * 11.1568e2;
2903 else if (x < proba+probhe+probt+probd+probp+probn+problamb0+probg) {
2914 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;
2918 else if (x < proba+probhe+probt+probd+probp+probn+problamb0+probg+probimf) {
2925 imf(af,zf,tcn,ee,&zimf,&aimf,&bimf,&sbimf,&timf,jprf);
2927 if(iloop>100)std::cout <<
"Problem in EVAPORA: IMF called > 100 times" << std::endl;
2928 if(zimf>=(zf-2.0))
goto dir1973;
2934 if(zimf==0.0 || aimf==0.0 || sbimf>ee)std::cout <<
"warning: Look in EVAPORA CALL IMF" << std::endl;
2945 tkeimf=
min(2.0*timf,ee-sbimf);
2948 if(tkeimf<=0.0)
goto dir1235;
2949 if(tkeimf>(ee-sbimf) && timf>0.5)
goto dir1235;
2951 tkeimf = tkeimf + bimf;
2955 epsiln = (sbimf-bimf) + tkeimf;
2984 if (ee <= 0.01)ee = 0.01;
2986 if(gammadecay==1 && ee<(epsiln+0.010)){
2987 epsiln = ee - 0.010;
2992 std::cout <<
"***WARNING epsilon<0***" << std::endl;
3000 if (ee <= 0.01)ee = 0.01;
3001 mtota = mtota + malpha;
3008 if(amoins==2 && zmoins==0){twon=1;amoins=1;}
else{ twon=0;}
3012 if(ff==0 && fimf==0){
3015 EV_TEMP[IEV_TAB_TEMP][0] = 0.;
3016 EV_TEMP[IEV_TAB_TEMP][1] = -2;
3017 EV_TEMP[IEV_TAB_TEMP][5] = 1.;
3019 EV_TEMP[IEV_TAB_TEMP][0] = zmoins;
3020 EV_TEMP[IEV_TAB_TEMP][1] = amoins;
3021 EV_TEMP[IEV_TAB_TEMP][5] = 0.;
3024 ctet1 = 2.0*rnd - 1.0;
3025 stet1 = std::sqrt(1.0 - std::pow(ctet1,2));
3027 phi1 = rnd*2.0*3.141592654;
3028 G4double xcv = stet1*std::cos(phi1);
3029 G4double ycv = stet1*std::sin(phi1);
3034 G4double ETOT_LP = std::sqrt(
pc*
pc + amoins*amoins * mu2);
3035 if(flamb0decay==1)ETOT_LP = std::sqrt(
pc*
pc + 1115.683*1115.683);
3036 EV_TEMP[IEV_TAB_TEMP][2] = c *
pc * xcv / ETOT_LP;
3037 EV_TEMP[IEV_TAB_TEMP][3] = c *
pc * ycv / ETOT_LP;
3038 EV_TEMP[IEV_TAB_TEMP][4] = c *
pc * zcv / ETOT_LP;
3041 EV_TEMP[IEV_TAB_TEMP][2] =
pc * xcv;
3042 EV_TEMP[IEV_TAB_TEMP][3] =
pc * ycv;
3043 EV_TEMP[IEV_TAB_TEMP][4] =
pc * zcv;
3045 G4double VXOUT=0.,VYOUT=0.,VZOUT=0.;
3047 EV_TEMP[IEV_TAB_TEMP][2],EV_TEMP[IEV_TAB_TEMP][3],
3048 EV_TEMP[IEV_TAB_TEMP][4],
3049 &VXOUT,&VYOUT,&VZOUT);
3050 EV_TEMP[IEV_TAB_TEMP][2] = VXOUT;
3051 EV_TEMP[IEV_TAB_TEMP][3] = VYOUT;
3052 EV_TEMP[IEV_TAB_TEMP][4] = VZOUT;
3055 G4double v2 = std::pow(EV_TEMP[IEV_TAB_TEMP][2],2.) +
3056 std::pow(EV_TEMP[IEV_TAB_TEMP][3],2.) +
3057 std::pow(EV_TEMP[IEV_TAB_TEMP][4],2.);
3058 G4double gamma = 1.0/std::sqrt(1.0 - v2 / (c*c));
3059 G4double etot_lp = amoins*mu * gamma;
3060 pxeva = pxeva - EV_TEMP[IEV_TAB_TEMP][2] * etot_lp / c;
3061 pyeva = pyeva - EV_TEMP[IEV_TAB_TEMP][3] * etot_lp / c;
3062 pleva = pleva - EV_TEMP[IEV_TAB_TEMP][4] * etot_lp / c;
3065 pxeva = pxeva - EV_TEMP[IEV_TAB_TEMP][2];
3066 pyeva = pyeva - EV_TEMP[IEV_TAB_TEMP][3];
3067 pleva = pleva - EV_TEMP[IEV_TAB_TEMP][4];
3069 G4double pteva = std::sqrt(pxeva*pxeva + pyeva*pyeva);
3071 G4double etot = std::sqrt ( pleva*pleva + pteva*pteva + af*af * mu2 );
3072 vxeva = c * pxeva / etot;
3073 vyeva = c * pyeva / etot;
3074 vleva = c * pleva / etot;
3075 IEV_TAB_TEMP = IEV_TAB_TEMP + 1;
3078 if(twon==1){
goto secondneutron;}
3081 if (zf < 3. || (ff == 1) || (fgamma == 1) || (fimf==1)) {
3093 (*tkeimf_par) = tkeimf;
3094 (*mtota_par) = mtota;
3095 (*vleva_par) = vleva;
3096 (*vxeva_par) = vxeva;
3097 (*vyeva_par) = vyeva;
3100 (*inttype_par) = inttype;
3101 (*iev_tab_temp_par)= IEV_TAB_TEMP;
3103 (*NbLam0_par) = NbLam0;
3107void 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)
3118 G4double problamb0 = (*problamb0_par);
3271 G4int choice_fisspart = 0;
3380 mglw(a-1.0,zprf,&ma1z);
3381 mglw(a-1.0,zprf-1.0,&ma1z1);
3382 mglw(a-2.0,zprf-1.0,&ma2z1);
3383 mglw(a-3.0,zprf-1.0,&ma3z1);
3384 mglw(a-3.0,zprf-2.0,&ma3z2);
3385 mglw(a-4.0,zprf-2.0,&ma4z2);
3388 if((a-1.)==3.0 && (zprf-1.0)==2.0) ma1z1=-7.7181660;
3389 if((a-1.)==4.0 && (zprf-1.0)==2.0) ma1z1=-28.295992;
3394 sd = ma2z1 - maz - 2.2246;
3395 st = ma3z1 - maz - 8.481977;
3396 she = ma3z2 - maz - 7.7181660;
3397 sa = ma4z2 - maz - 28.295992;
3427 if (zprf <= 1.0e0 || a <= 1.0e0 || (a-zprf) < 0.0) {
3437 if (zprf <= 1.0e0 || a <= 2.0e0 || (a-zprf) < 1.0) {
3447 if (zprf <= 1.0e0 || a <= 3.0e0 || (a-zprf) < 2.0) {
3457 if (a-4.0<=0.0 || zprf<=2.0 || (a-zprf)<2.0) {
3467 if (a-3.0 <= 0.0 || zprf<=2.0 || (a-zprf)<1.0) {
3481 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);
3492 barfit(k,k+j,il,&sbfis,&segs,&selmax);
3498 ef = ef*(4.5114-2.2687*(a-zprf)/zprf);
3500 ef = ef*(3.3931-1.5338*(a-zprf)/zprf);
3504 if((a-zprf)/zprf>1.52)ef=ef*(1.1222-0.10886*(a-zprf)/zprf)-0.1;
3506 if(k>=94&&k<=98&&j<158){
3508 if(
mod(j,2)==0&&
mod(k,2)==0){
3509 if(k>=94){ef = ef-(11.54108*(a-zprf)/zprf-18.074);}
3512 if(
mod(j,2)==1&&
mod(k,2)==1){
3513 if(k>=95){ef = ef-(14.567*(a-zprf)/zprf-23.266);}
3516 if(
mod(j,2)==0&&
mod(k,2)==1){
3517 if(j>=144){ef = ef-(13.662*(a-zprf)/zprf-21.656);}
3520 if(
mod(j,2)==1&&
mod(k,2)==0){
3521 if(j>=144){ef = ef-(13.662*(a-zprf)/zprf-21.656);}
3532 if (ef < 0.0)ef = 0.0;
3538 ef = ef + 0.51*(1115.-938.+sn-slamb0)/std::pow(a,2./3.);
3571 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);
3572 erot = jprf * jprf * 197.328 * 197.328 /(2. * iinert);
3575 bsbkbc(a,zprf,&bscn,&bkcn,&bccn);
3578 densniv(a,zprf,ee,0.0,&densg,bshell,bscn,bkcn,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprf,0,&qrcn);
3630 lorb(a,a-1.,jprf,ee-sn,&dlout,&sdlout);
3632 if(IDjprf==1) djprf = 0.0;
3633 jprfn = jprf + djprf;
3634 jprfn =
dint(std::abs(jprfn));
3639 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);
3640 erotn = jprfn * jprfn * 197.328 * 197.328 /(2. * iinert);
3641 bsbkbc(a-1.,zprf,&bs,&bk,&bc);
3644 densniv(a-1.0,zprf,ee,sn,&densn,bshell, bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprfn,0,&qr);
3654 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ_NEUT CALLED MORE THAN 100 TIMES" << std::endl;
3663 if(ecn<=0.0)
goto dir1234;
3683 lorb(a,a-1.,jprf,ee-sbp,&dlout,&sdlout);
3685 if(IDjprf==1) djprf = 0.0;
3686 jprfp = jprf + djprf;
3687 jprfp =
dint(std::abs(jprfp));
3692 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);
3693 erotp = jprfp * jprfp * 197.328 * 197.328 /(2. * iinert);
3695 bsbkbc(a-1.,zprf-1.,&bs,&bk,&bc);
3698 densniv(a-1.0,zprf-1.0,ee,sbp,&densp,bshell,bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprfp,0,&qr);
3708 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3717 if(ecp<=0.0)
goto dir1235;
3720 ecp = 2.0 * pt +
bp;
3734 if ((in >= 2) && (iz >= 2)) {
3738 lorb(a,a-2.,jprf,ee-sbd,&dlout,&sdlout);
3740 if(IDjprf==1) djprf = 0.0;
3741 jprfd = jprf + djprf;
3742 jprfd =
dint(std::abs(jprfd));
3747 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);
3748 erotd = jprfd * jprfd * 197.328 * 197.328 /(2. * iinert);
3750 bsbkbc(a-2.,zprf-1.,&bs,&bk,&bc);
3753 densniv(a-2.0,zprf-1.0e0,ee,sbd,&densd,bshell,bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprfd,0,&qr);
3764 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3773 if(ecd<=0.0)
goto dir1236;
3776 ecd = 2.0 * dt + bd;
3790 if ((in >= 3) && (iz >= 2)) {
3794 lorb(a,a-3.,jprf,ee-sbt,&dlout,&sdlout);
3796 if(IDjprf==1) djprf = 0.0;
3797 jprft = jprf + djprf;
3798 jprft =
dint(std::abs(jprft));
3803 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);
3804 erott = jprft * jprft * 197.328 * 197.328 /(2. * iinert);
3806 bsbkbc(a-3.,zprf-1.,&bs,&bk,&bc);
3809 densniv(a-3.0,zprf-1.0,ee,sbt,&denst,bshell,bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprft,0,&qr);
3820 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3829 if(ect<=0.0)
goto dir1237;
3832 ect = 2.0 * tt + bt;
3846 if ((in >= 3) && (iz >= 3)) {
3850 lorb(a,a-4.,jprf,ee-sba,&dlout,&sdlout);
3852 if(IDjprf==1) djprf = 0.0;
3853 jprfa = jprf + djprf;
3854 jprfa =
dint(std::abs(jprfa));
3859 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);
3860 erota = jprfa * jprfa * 197.328 * 197.328 /(2. * iinert);
3862 bsbkbc(a-4.,zprf-2.,&bs,&bk,&bc);
3865 densniv(a-4.0,zprf-2.0,ee,sba,&densa,bshell,bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprfa,0,&qr);
3876 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3885 if(eca<=0.0)
goto dir1238;
3888 eca = 2.0 * at + ba;
3902 if ((in >= 2) && (iz >= 3)) {
3906 lorb(a,a-3.,jprf,ee-sbhe,&dlout,&sdlout);
3908 if(IDjprf==1) djprf = 0.0;
3909 jprfhe = jprf + djprf;
3910 jprfhe =
dint(std::abs(jprfhe));
3915 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);
3916 erothe = jprfhe * jprfhe * 197.328 * 197.328 /(2. * iinert);
3918 bsbkbc(a-3.,zprf-2.,&bs,&bk,&bc);
3921 densniv(a-3.0,zprf-2.0,ee,sbhe,&denshe,bshell,bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprfhe,0,&qr);
3932 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ CALLED MORE THAN 100 TIMES" << std::endl;
3941 if(eche<=0.0)
goto dir1239;
3944 eche = 2.0 * het + bhe;
3960 if (in >= 2 && NbLam0>0) {
3964 lorb(a,a-1.,jprf,ee-slamb0,&dlout,&sdlout);
3966 if(IDjprf==1) djprf = 0.0;
3967 jprflamb0 = jprf + djprf;
3968 jprflamb0 =
dint(std::abs(jprflamb0));
3973 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);
3974 erotlamb0 = jprflamb0 * jprflamb0 * 197.328 * 197.328 /(2. * iinert);
3975 bsbkbc(a-1.,zprf,&bs,&bk,&bc);
3978 densniv(a-1.0,zprf,ee,slamb0,&denslamb0,bshell, bs,bk,&temp,
fiss->
optshp,
fiss->
optcol,defbet,&ecor,jprflamb0,0,&qr);
3988 if(IS>100){std::cout <<
"WARNING: FVMAXHAZ_NEUT CALLED MORE THAN 100 TIMES" << std::endl;
3991 if(eclamb0>(ee-slamb0)){
3992 if((ee-slamb0)<rlamb0t)
3993 eclamb0 = ee-slamb0;
3997 if(eclamb0<=0.0)
goto dir1240;
3999 eclamb0 = 2.0 * lamb0t;
4027 gn =
width(a,zprf,1.0,0.0,nt,0.0,sn,ee-erotn)* densn/densg;
4032 gp =
width(a,zprf,1.0,1.0,pt,
bp,sbp,ee-erotp)*densp/densg*
pen(a, 1.0, omegap, pt);
4037 gd =
width(a,zprf,2.0,1.0,dt,bd,sbd,ee-erotd)*densd/densg*
pen(a, 2.0, omegad, dt);
4042 gt =
width(a,zprf,3.0,1.0,tt,bt,sbt,ee-erott)*denst/densg*
pen(a, 3.0, omegat, tt);
4047 ghe =
width(a,zprf,3.0,2.0,het,bhe,sbhe,ee-erothe) * denshe/densg*
pen(a, 3.0, omegahe, het);
4052 ga =
width(a,zprf,4.0,2.0,at,ba,sba,ee-erota) * densa/densg*
pen(a, 4.0, omegaa, at);
4057 glamb0 =
width(a,zprf,1.0,-2.0,lamb0t,0.0,slamb0,ee-erotlamb0)* denslamb0/densg;
4065 G4int izcn=0,incn=0,inmin=0,inmax=0,inmi=0,inma=0;
4068 if(fimf_allowed==0 || zprf<=5.0 || a<=7.0){
4086 inmin =
max(inmin,incn-inma);
4087 inmax =
min(inmax,incn-inmi);
4089 inmax =
max(inmax,inmin);
4091 for(
G4int iaimf=izimf+inmin;iaimf<=izimf+inmax;iaimf++){
4093 if(aimf>=a || zimf>=zprf){
4105 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));
4107 erot = jprf * jprf * 197.328 * 197.328 /(2.0 * iinert);
4110 if(densg==0.0 || ee < (sbimf + erot)){
4116 densniv(a,zprf,ee,sbimf,&densimf,0.0,bsimf,1.0,&timf,0,0,defbetimf,&ecor,jprf,2,&qr);
4118 imfarg = (sbimf+erotcn-erot)/timf;
4119 if(imfarg > 200.0) imfarg = 200.0;
4129 width_imf =
width(a,zprf,aimf,zimf,timf,bimf,sbimf,ee-erot)*std::exp(-imfarg)*qr/qrcn;
4132 gimf3 = gimf3 + width_imf;
4146 inmin =
max(inmin,incn-inma);
4147 inmax =
min(inmax,incn-inmi);
4149 inmax =
max(inmax,inmin);
4151 for(
G4int iaimf=izimf+inmin;iaimf<=izimf+inmax;iaimf++){
4153 if(aimf>=a || zimf>=zprf){
4165 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));
4167 erot = jprf * jprf * 197.328 * 197.328 /(2.0 * iinert);
4170 if(densg==0.0 || ee < (sbimf + erot)){
4176 densniv(a,zprf,ee,sbimf,&densimf,0.0,bsimf,1.0,&timf,0,0,defbetimf,&ecor,jprf,2,&qr);
4178 imfarg = (sbimf+erotcn-erot)/timf;
4179 if(imfarg > 200.0) imfarg = 200.0;
4188 width_imf =
width(a,zprf,aimf,zimf,timf,bimf,sbimf,ee-erot)*std::exp(-imfarg)*qr/qrcn;
4191 gimf5 = gimf5 + width_imf;
4196 if(gimf3<=0.0 || gimf5<=0.0){
4202 b_imf = (std::log10(gimf3) - std::log10(gimf5))/(std::log10(3.0)-std::log10(5.0));
4204 if(b_imf >= -1.01) b_imf = -1.01;
4205 if(b_imf <= -100.0) {
4212 a_imf = gimf3 / std::pow(3.0,b_imf);
4213 gimf = a_imf * ( std::pow(zprf,b_imf+1.0) - std::pow(3.0,b_imf+1.0)) /(b_imf + 1.0);
4217 if(gimf < 1.e-10) gimf = 0.0;
4224 pa = (
ald->
av)*a + (
ald->
as)*std::pow(a,2./3.) + (
ald->
ak)*std::pow(a,1./3.);
4225 gamma = 2.5 * pa * std::pow(a,-4./3.);
4231 gtemp = 17.60/(std::pow(a,0.699) * std::sqrt(gfactor));
4234 gg = 0.624e-9*std::pow(a,1.6)*std::pow(gtemp,5.);
4244 gsum = ga + ghe + gd + gt + gp + gn + gimf + gg + glamb0;
4257 if(
fiss->
ifis==0 || (zprf*zprf/a<=22.74 && zprf<60.)){
4265 fission_width(zprf,a,ee,bssp,bksp,ef,y,&gf,&temp,jprf,0,1,
fiss->
optcol,
fiss->
optshp,densg);
4285 gtotal = ga + ghe + gp + gd + gt + gn + gg +gimf + gf + glamb0;
4304 probhe = ghe/gtotal;
4307 probimf = gimf/gtotal;
4308 problamb0 = glamb0/gtotal;
4321 fomega_sp(a,y,&mfcd,&omegasp,&homegasp);
4325 fomega_gs(a,zprf,&k1,&omegags,&homegags);
4338 part_fiss(
fiss->
bet,gsum,gf,y,tauc,ts1,tsum, &choice_fisspart,zprf,a,ft,&t_lapse,&gff);
4342 tsum = tsum + t_lapse;
4345 if(choice_fisspart==2){
4361 gtotal=ga + ghe + gp + gd + gt + gn + gimf + gg + glamb0;
4380 probhe = ghe/gtotal;
4383 probimf = gimf/gtotal;
4384 problamb0 = glamb0/gtotal;
4390 gtotal = ga + ghe + gp + gd + gt + gn + gg + gimf + glamb0;
4408 probhe = ghe/gtotal;
4411 probimf = gimf/gtotal;
4412 problamb0 = glamb0/gtotal;
4416 ptotl = probp+probd+probt+probn+probhe+proba+probg+probimf+probf+problamb0;
4422 (*probp_par) = probp;
4423 (*probd_par) = probd;
4424 (*probt_par) = probt;
4425 (*probn_par) = probn;
4426 (*probhe_par) = probhe;
4427 (*proba_par) = proba;
4428 (*probg_par) = probg;
4429 (*probimf_par) = probimf;
4430 (*problamb0_par) = problamb0;
4431 (*probf_par) = probf;
4432 (*ptotl_par) = ptotl;
4439 (*slamb0_par) = slamb0;
4452 (*eclamb0_par) = eclamb0;
4460 (*jprfn_par) = jprfn;
4461 (*jprfp_par) = jprfp;
4462 (*jprfd_par) = jprfd;
4463 (*jprft_par) = jprft;
4464 (*jprfhe_par) = jprfhe;
4465 (*jprfa_par) = jprfa;
4466 (*jprflamb0_par) = jprflamb0;
4471void 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)
4566 G4double pi6 = std::pow(3.1415926535,2) / 6.0;
4578 if(afp<=20) BSHELLCT = 0.0;
4611 pa = (
ald->
av)*a + (
ald->
as)*std::pow(a,(2.e0/3.e0)) + (
ald->
ak)*std::pow(a,(1.e0/3.e0));
4613 pa = (
ald->
av)*a + (
ald->
as)*bsin*std::pow(a,(2.e0/3.e0)) + (
ald->
ak)*bkin*std::pow(a,(1.e0/3.e0));
4615 gamma = 2.5 * pa * std::pow(a,-4.0/3.0);
4620 if(ifis==0&&bs!=1.0){
4623 if(ponq>700.0) ponq = 700.0;
4624 bs = 1.0/(1.0+std::exp(-ponq)) + 1.0/(1.0+std::exp(ponq)) * bsin;
4625 bk = 1.0/(1.0+std::exp(-ponq)) + 1.0/(1.0+std::exp(ponq)) * bkin;
4630 pa = (
ald->
av)*a + (
ald->
as)*std::pow(a,(2.e0/3.e0)) + (
ald->
ak)*std::pow(a,(1.e0/3.e0));
4633 pa = (
ald->
av)*a + (
ald->
as)*bs*std::pow(a,(2.e0/3.e0)) + (
ald->
ak)*bk*std::pow(a,(1.e0/3.e0));
4636 gamma = 2.5 * pa * std::pow(a,-4.0/3.0);
4642 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));
4662 deltpp = -0.25e0* std::pow((delta0/std::sqrt(a)),2) * pa /pi6 + 22.34e0*std::pow(a,-0.464)-0.235;
4666 e=e-(0.285+11.17*std::pow(a,-0.464)-0.390-0.00058*a);
4670 e=e-(22.34*std::pow(a,-0.464)-0.235);
4694 ponfe = -2.5*pa*e*std::pow(a,(-4.0/3.0));
4696 if (ponfe < -700.0) {
4699 fe = 1.0 - std::exp(ponfe);
4702 he = 1.0 - std::pow((1.0 - e/ecr),2);
4709 fecor = e + deltau*
fe + deltpp*he;
4716 y1 = std::sqrt(pa*fecor);
4717 for(
G4int j = 0; j < 5; j++) {
4718 y2 = pa*fecor*(1.e0-std::exp(-y1));
4723 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;
4726 y11 = std::sqrt(pa*ecor1);
4727 for(
G4int j = 0; j < 7; j++) {
4728 y21 = pa*ecor1*(1.0-std::exp(-y11));
4729 y11 = std::sqrt(y21);
4733 fdens = fdens*std::pow((y01/y0),1.5);
4734 ftemp = ftemp*std::pow((y01/y0),1.5);
4738 ponniv = 2.0*std::sqrt(pa*fecor);
4739 if (ponniv > 700.0) {
4743 fdens = 0.1477045 * std::exp(ponniv)/(std::pow(pa,0.25)*std::pow(fecor,1.25));
4744 ftemp = std::sqrt(fecor/pa);
4750 if(IOPTCT==0)
goto densniv100;
4751 tempct = 17.60/( std::pow(a,0.699) * std::sqrt(1.+gamma*BSHELLCT));
4762 if (IPARITE == 1) { e0 = 0.285+11.17*std::pow(a,-0.464) - 0.390-0.00058*a;}
4764 if (IPARITE == 2) { e0 = 22.34*std::pow(a,-0.464)-0.235;}
4766 if (IPARITE == 0){ e0 = 0.0;}
4768 ponniv = (ein-e0)/tempct;
4769 if(ifis!=1) ponniv =
max(0.0,(ein-e0)/tempct);
4770 if(ponniv>700.0){ ponniv = 700.0;}
4771 densct = std::exp(ponniv)/tempct*std::exp(0.079*BSHELLCT/tempct);
4775 if(elim>=ecr&&densfm<=densct){
4783 if(elim>=ecr&&tfm>=tempct){
4791 ponniv = (ein)/tempct;
4792 if(ponniv>700.0){ ponniv = 700.0;}
4793 densct = std::exp(ponniv)/tempct;
4795 if(ein>=ecr && densfm<=densct){
4805 if(ein>=ecr && tfm>=tempct){
4820 ftemp = 17.60/( std::pow(a,0.699) * std::sqrt(1.0+gamma*BSHELLCT));
4833 fnorm = std::pow(1.16,2)*931.49*1.e-2/(9.0* std::pow(6.582122,2));
4835 if(ifis==0 || ifis==2){
4841 fp_per = 0.4*std::pow(a,5.0/3.0)*fnorm*(1.0+0.50*defbet*std::sqrt(5.0/(4.0*
pi)));
4842 fp_par = 0.40*std::pow(a,5.0/3.0)*fnorm*(1.0-defbet*std::sqrt(5.0/(4.0*
pi)));
4851 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));
4853 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));
4860 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);
4861 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);
4865 if(fp_par<0.0)fp_par=0.0;
4866 if(fp_per<0.0)fp_per=0.0;
4868 sig_per = std::sqrt(fp_per * ftemp);
4869 sig_par = std::sqrt(fp_par * ftemp);
4871 sigma2 = sig_per*sig_per + sig_par*sig_par;
4872 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));
4873 erot = jprf*jprf/(2.0*std::sqrt(fp_par*fp_par+fp_per*fp_per));
4877 qrot(z,a,defbet,sig_per,fecor-erot,&fqr);
4883 fdens = fdens * fqr *jfact;
4885 if(fdens<1e-300)fdens=0.0;
4921 G4double ponq = 0.0, dn = 0.0,
n = 0.0, dz = 0.0;
4922 G4int distn,distz,ndist, zdist;
4923 G4int nmn[8]= {2, 8, 14, 20, 28, 50, 82, 126};
4924 G4int nmz[8]= {2, 8, 14, 20, 28, 50, 82, 126};
4928 if(std::abs(bet)<=0.15){
4939 for(
G4int i =0;i<8;i++){
4940 ndist = std::fabs(
idnint(
n) - nmn[i]);
4941 if(ndist < distn) distn = ndist;
4942 zdist = std::fabs(
idnint(z) - nmz[i]);
4943 if(zdist < distz) distz = zdist;
4949 bet = 0.022 + 0.003*dn + 0.002*dz;
4951 sig = 75.0*std::pow(bet,2.) * sig;
4955 ponq = (u - ucr)/dcr;
4963 (*qr) = 1.0/(1.0 + std::exp(ponq)) * (sig - 1.0) + 1.0;
4984 for(
G4int i = 2; i <
n; i++) {
5010 if(ia==0)
return eflmacResult;
5013 G4double z = 0.0,
n = 0.0, a = 0.0, av = 0.0, as = 0.0;
5014 G4double a0 = 0.0, c1 = 0.0, c4 = 0.0, b1 = 0.0, b3 = 0.0;
5015 G4double ff = 0.0, ca = 0.0, w = 0.0, efl = 0.0;
5016 G4double r0 = 0.0, kf = 0.0, ks = 0.0;
5017 G4double kv = 0.0, rp = 0.0, ay = 0.0, aden = 0.0, x0 = 0.0, y0 = 0.0;
5018 G4double esq = 0.0, ael = 0.0, i = 0.0, e0 = 0.0;
5068 if(flag==1){
goto eflmac311;}
5078 c1 = 3.0/5.0*esq/r0;
5079 c4 = 5.0/4.0*std::pow((3.0/(2.0*
pi)),(2.0/3.0)) * c1;
5080 kf = std::pow((9.0*
pi*z/(4.0*a)),(1.0/3.0))/r0;
5082 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));
5085 x0 = r0 * std::pow(a,(1.0/3.0)) / ay;
5086 y0 = r0 * std::pow(a,(1.0/3.0)) / aden;
5088 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);
5090 b3 = 1.0 - 5.0/std::pow(y0,2) * (1.0 - 15.0/(8.0*y0) + 21.0/(8.0 * std::pow(y0,3))
5091 - 3.0/4.0 * (1.0 + 9.0/(2.0*y0) + 7.0/std::pow(y0,2)
5092 + 7.0/(2.0 * std::pow(y0,3))) * std::exp(-2.0*y0));
5096 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
5097 + 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))
5098 + ff*std::pow(z,2)/a -ca*(
n-z) - ael * std::pow(z,(2.39e0));
5100 efl = efl + w*std::abs(i);
5105 if (in==iz && (
mod(in,2) == 1) && (
mod(iz,2) == 1) && in>0.) {
5121 e0 = 0.285+11.17*std::pow(a,-0.464) -0.390-0.00058*(a);
5127 e0 = 22.34*std::pow(a,-0.464)-0.235;
5139 return eflmacResult;
5164 (*del) = -12.0/std::sqrt(a);
5167 (*del) = 12.0/std::sqrt(a);
5179 G4double n1 = 0.0, n2 = 0.0, n3 = 0.0;
5215 if (bet/(std::sqrt(2.0)*10.0*(homega/6.582122)) <= 1.0) {
5216 tauResult = std::log(10.0*ef/t)/(bet*1.0e21);
5219 tauResult = std::log(10.0*ef/t)/ (2.0*std::pow((10.0*homega/6.582122),2))*(bet*1.0e-21);
5231 G4double rel = bet/(20.0*homega/6.582122);
5232 G4double cramResult = std::sqrt(1.0 + std::pow(rel,2)) - rel;
5235 if (cramResult > 1.0) {
5257 const G4int bsbkSize = 54;
5259 G4double bk[bsbkSize] = {0.0, 1.00000,1.00087,1.00352,1.00799,1.01433,1.02265,1.03306,
5260 1.04576,1.06099,1.07910,1.10056,1.12603,1.15651,1.19348,
5261 1.23915,1.29590,1.35951,1.41013,1.44103,1.46026,1.47339,
5262 1.48308,1.49068,1.49692,1.50226,1.50694,1.51114,1.51502,
5263 1.51864,1.52208,1.52539,1.52861,1.53177,1.53490,1.53803,
5264 1.54117,1.54473,1.54762,1.55096,1.55440,1.55798,1.56173,
5265 1.56567,1.56980,1.57413,1.57860,1.58301,1.58688,1.58688,
5266 1.58688,1.58740,1.58740, 0.0};
5268 G4double bs[bsbkSize] = {0.0, 1.00000,1.00086,1.00338,1.00750,1.01319,
5269 1.02044,1.02927,1.03974,
5270 1.05195,1.06604,1.08224,1.10085,1.12229,1.14717,1.17623,1.20963,
5271 1.24296,1.26532,1.27619,1.28126,1.28362,1.28458,1.28477,1.28450,
5272 1.28394,1.28320,1.28235,1.28141,1.28042,1.27941,1.27837,1.27732,
5273 1.27627,1.27522,1.27418,1.27314,1.27210,1.27108,1.27006,1.26906,
5274 1.26806,1.26707,1.26610,1.26514,1.26418,1.26325,1.26233,1.26147,
5275 1.26147,1.26147,1.25992,1.25992, 0.0};
5277 i =
idint(y/(2.0e-02)) + 1;
5279 if((i + 1) >= bsbkSize) {
5287 bipolResult = bs[i] + (bs[i+1] - bs[i])/2.0e-02 * (y - 2.0e-02*(i - 1));
5290 bipolResult = bk[i] + (bk[i+1] - bk[i])/2.0e-02 * (y - 2.0e-02*(i - 1));
5305 ES0 = 20.760*std::pow(AF,2.0/3.0);
5308 MR02 = std::pow(AF,5.0/3.0)*1.0340*0.010*1.175*1.175;
5310 (*MFCD) = MR02 * 3.0/10.0*(1.0+3.0*
Y);
5312 OMEGA = std::sqrt(ES0/MR02)*std::sqrt(8.0/3.0*
Y*(1.0+304.0*
Y/255.0));
5314 HOMEGA = 6.58122*OMEGA/10.0;
5329 G4double OMEGA,HOMEGA,MR02,MINERT,
C,fk1;
5331 MR02 = std::pow(AF,5.0/3.0)*1.0340*0.01*1.175*1.175;
5332 MINERT = 3.*MR02/10.0;
5333 C = 17.9439*(1.-1.7826*std::pow((AF-2.0*ZF)/AF,2));
5334 fk1 = 0.4*
C*std::pow(AF,2.0/3.0)-0.1464*std::pow(ZF,2)/std::pow(AF,1./3.);
5335 OMEGA = std::sqrt(fk1/MINERT);
5336 HOMEGA = 6.58122*OMEGA/10.0;
5371 BARR = 1.345 *
Z1 * Z2 / RMAX;
5373 OMEGA = 4.5 / 197.3287;
5457 for(
G4int init_i = 0; init_i < 7; init_i++) {
5461 for(
G4int init_i = 0; init_i < 10; init_i++) {
5466 G4double amax2 = 0.0, aa = 0.0, zz = 0.0, bfis = 0.0;
5467 G4double bfis0 = 0.0, ell = 0.0, el = 0.0, egs = 0.0, el80 = 0.0, el20 = 0.0;
5468 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;
5469 G4double aj = 0.0, ak = 0.0, a1 = 0.0, a2 = 0.0;
5471 G4int i = 0, j = 0, k = 0,
m = 0;
5474 G4double emncof[4][5] = {{-9.01100e+2,-1.40818e+3, 2.77000e+3,-7.06695e+2, 8.89867e+2},
5475 {1.35355e+4,-2.03847e+4, 1.09384e+4,-4.86297e+3,-6.18603e+2},
5476 {-3.26367e+3, 1.62447e+3, 1.36856e+3, 1.31731e+3, 1.53372e+2},
5477 {7.48863e+3,-1.21581e+4, 5.50281e+3,-1.33630e+3, 5.05367e-2}};
5479 G4double elmcof[4][5] = {{1.84542e+3,-5.64002e+3, 5.66730e+3,-3.15150e+3, 9.54160e+2},
5480 {-2.24577e+3, 8.56133e+3,-9.67348e+3, 5.81744e+3,-1.86997e+3},
5481 {2.79772e+3,-8.73073e+3, 9.19706e+3,-4.91900e+3, 1.37283e+3},
5482 {-3.01866e+1, 1.41161e+3,-2.85919e+3, 2.13016e+3,-6.49072e+2}};
5484 G4double emxcof[4][6] = {{9.43596e4,-2.241997e5,2.223237e5,-1.324408e5,4.68922e4,-8.83568e3},
5485 {-1.655827e5,4.062365e5,-4.236128e5,2.66837e5,-9.93242e4,1.90644e4},
5486 {1.705447e5,-4.032e5,3.970312e5,-2.313704e5,7.81147e4,-1.322775e4},
5487 {-9.274555e4,2.278093e5,-2.422225e5,1.55431e5,-5.78742e4,9.97505e3}};
5489 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},
5490 {-1.13269453e+6, 2.97764590e+6,-4.54326326e+6, 3.00464870e+6, -1.44989274e+6,-1.02026610e+5, 6.27959815e+4},
5491 {1.37543304e+6,-3.65808988e+6, 5.47798999e+6,-3.78109283e+6, 1.84131765e+6, 1.53669695e+4,-6.96817834e+4},
5492 {-8.56559835e+5, 2.48872266e+6,-4.07349128e+6, 3.12835899e+6, -1.62394090e+6, 1.19797378e+5, 4.25737058e+4},
5493 {3.28723311e+5,-1.09892175e+6, 2.03997269e+6,-1.77185718e+6, 9.96051545e+5,-1.53305699e+5,-1.12982954e+4},
5494 {4.15850238e+4, 7.29653408e+4,-4.93776346e+5, 6.01254680e+5, -4.01308292e+5, 9.65968391e+4,-3.49596027e+3},
5495 {-1.82751044e+5, 3.91386300e+5,-3.03639248e+5, 1.15782417e+5, -4.24399280e+3,-6.11477247e+3, 3.66982647e+2}};
5497 const G4int sizex = 5;
5498 const G4int sizey = 6;
5499 const G4int sizez = 4;
5501 G4double egscof[sizey][sizey][sizez];
5503 G4double egs1[sizey][sizex] = {{1.927813e5, 7.666859e5, 6.628436e5, 1.586504e5,-7.786476e3},
5504 {-4.499687e5,-1.784644e6,-1.546968e6,-4.020658e5,-3.929522e3},
5505 {4.667741e5, 1.849838e6, 1.641313e6, 5.229787e5, 5.928137e4},
5506 {-3.017927e5,-1.206483e6,-1.124685e6,-4.478641e5,-8.682323e4},
5507 {1.226517e5, 5.015667e5, 5.032605e5, 2.404477e5, 5.603301e4},
5508 {-1.752824e4,-7.411621e4,-7.989019e4,-4.175486e4,-1.024194e4}};
5510 G4double egs2[sizey][sizex] = {{-6.459162e5,-2.903581e6,-3.048551e6,-1.004411e6,-6.558220e4},
5511 {1.469853e6, 6.564615e6, 6.843078e6, 2.280839e6, 1.802023e5},
5512 {-1.435116e6,-6.322470e6,-6.531834e6,-2.298744e6,-2.639612e5},
5513 {8.665296e5, 3.769159e6, 3.899685e6, 1.520520e6, 2.498728e5},
5514 {-3.302885e5,-1.429313e6,-1.512075e6,-6.744828e5,-1.398771e5},
5515 {4.958167e4, 2.178202e5, 2.400617e5, 1.167815e5, 2.663901e4}};
5517 G4double egs3[sizey][sizex] = {{3.117030e5, 1.195474e6, 9.036289e5, 6.876190e4,-6.814556e4},
5518 {-7.394913e5,-2.826468e6,-2.152757e6,-2.459553e5, 1.101414e5},
5519 {7.918994e5, 3.030439e6, 2.412611e6, 5.228065e5, 8.542465e3},
5520 {-5.421004e5,-2.102672e6,-1.813959e6,-6.251700e5,-1.184348e5},
5521 {2.370771e5, 9.459043e5, 9.026235e5, 4.116799e5, 1.001348e5},
5522 {-4.227664e4,-1.738756e5,-1.795906e5,-9.292141e4,-2.397528e4}};
5524 G4double egs4[sizey][sizex] = {{-1.072763e5,-5.973532e5,-6.151814e5, 7.371898e4, 1.255490e5},
5525 {2.298769e5, 1.265001e6, 1.252798e6,-2.306276e5,-2.845824e5},
5526 {-2.093664e5,-1.100874e6,-1.009313e6, 2.705945e5, 2.506562e5},
5527 {1.274613e5, 6.190307e5, 5.262822e5,-1.336039e5,-1.115865e5},
5528 {-5.715764e4,-2.560989e5,-2.228781e5,-3.222789e3, 1.575670e4},
5529 {1.189447e4, 5.161815e4, 4.870290e4, 1.266808e4, 2.069603e3}};
5531 for(i = 0; i < sizey; i++) {
5532 for(j = 0; j < sizex; j++) {
5533 egscof[i][j][0] = egs1[i][j];
5534 egscof[i][j][1] = egs2[i][j];
5535 egscof[i][j][2] = egs3[i][j];
5536 egscof[i][j][3] = egs4[i][j];
5541 if (iz < 19 || iz > 111) {
5545 if(iz > 102 && il > 0) {
5552 amin= 1.2e0*z + 0.01e0*z*z;
5553 amax= 5.8e0*z - 0.024e0*z*z;
5555 if(a < amin || a >
amax) {
5567 for(i = 0; i < 7; i++) {
5568 for(j = 0; j < 7; j++) {
5569 bfis0=bfis0+elzcof[j][i]*pz[i]*pa[j];
5581 amin2 = 1.4e0*z + 0.009e0*z*z;
5582 amax2 = 20.e0 + 3.0e0*z;
5584 if((a < amin2-5.e0 || a > amax2+10.e0) && il > 0) {
5594 for(i = 0; i < 4; i++) {
5595 for(j = 0; j < 5; j++) {
5596 el80 = el80 + elmcof[i][j]*pz[j]*pa[i];
5597 el20 = el20 + emncof[i][j]*pz[j]*pa[i];
5608 for(i = 0; i < 4; i++) {
5609 for(j = 0; j < 6; j++) {
5610 elmax = elmax + emxcof[i][j]*pz[j]*pa[i];
5621 x = sel20/(*selmax);
5622 y = sel80/(*selmax);
5626 q = 0.2/(std::pow(sel20,2)*std::pow(sel80,2)*(sel20-sel80));
5627 qa = q*(4.0*std::pow(sel80,3) - std::pow(sel20,3));
5628 qb = -q*(4.0*std::pow(sel80,2) - std::pow(sel20,2));
5629 bfis = bfis*(1.0 + qa*std::pow(el,2) + qb*std::pow(el,3));
5633 aj = (-20.0*std::pow(x,5) + 25.e0*std::pow(x,4) - 4.0)*std::pow((y-1.0),2)*y*y;
5634 ak = (-20.0*std::pow(y,5) + 25.0*std::pow(y,4) - 1.0) * std::pow((x-1.0),2)*x*x;
5635 q = 0.2/(std::pow((y-x)*((1.0-x)*(1.0-y)*x*y),2));
5636 qa = q*(aj*y - ak*x);
5637 qb = -q*(aj*(2.0*y + 1.0) - ak*(2.0*x + 1.0));
5639 a1 = 4.0*std::pow(z,5) - 5.0*std::pow(z,4) + 1.0;
5640 a2 = qa*(2.e0*z + 1.e0);
5641 bfis=bfis*(a1 + (z - 1.e0)*(a2 + qb*z)*z*z*(z - 1.e0));
5648 if(el > (*selmax)) {
5654 if(el > (*selmax)) {
5658 for(k = 0; k < 4; k++) {
5659 for(l = 0; l < 6; l++) {
5660 for(
m = 0;
m < 5;
m++) {
5661 egs = egs + egscof[l][
m][k]*pz[l]*pa[k]*pl[2*
m];
5707 ferf=-
gammp(0.5,x*x);
5709 ferf=
gammp(0.5,x*x);;
5719 if(x<0.0 || a<=0.0)std::cout <<
"G4Abla::gammp = bad arguments in gammp" << std::endl;
5721 gser(&gamser,a,x,gln);
5724 gcf(&gammcf,a,x,gln);
5743 for(
G4int i=1;i<=itmax;i++){
5747 if(std::fabs(d)<fpmin)d=fpmin;
5749 if(std::fabs(c)<fpmin)c=fpmin;
5753 if(std::fabs(del-1.)<
eps)
goto dir1;
5755 std::cout <<
"a too large, ITMAX too small in gcf" << std::endl;
5757 fgammcf=std::exp(-x+a*std::log(x)-gln)*h;
5770 if(x<0.)std::cout <<
"G4Abla::gser = x < 0 in gser" << std::endl;
5781 if(std::fabs(del)<std::fabs(sum)*
eps)
goto dir1;
5783 std::cout <<
"a too large, ITMAX too small in gser" << std::endl;
5785 fgamser=sum*std::exp(-x+a*std::log(x)-gln);
5793 G4double cof[6]={76.18009172947146,-86.50532032941677,24.01409824083091,
5794-1.231739572450155,0.1208650973866179e-2,-0.5395239384953e-5};
5800 tmp=(x+0.5)*std::log(tmp)-tmp;
5801 ser=1.000000000190015;
5802 for(
G4int j=0;j<6;j++){
5807 return fgammln=tmp+std::log(stp*ser/x);
5815 return (E*std::exp(-E));
5821 return (1.0 - (E + 1.0) * std::exp(-E));
5837 const G4int pSize = 101;
5855 for(i = 1; i <= 99; i++) {
5859 if (std::fabs(
f(x) -
G4double(i)/100.0) < 1e-5) {
5884 x = (p[i] - p[i-1])*(y*100 - i) + p[i];
5901 if(ii <= 0 || jj < 0) {
5912 fpace2=fpace2/1000.;
5914 if(
pace->
dm[ii][jj] == 0.) {
5919 guet(&a, &z, &fpace2);
5920 fpace2=fpace2-ii*931.5;
5921 fpace2=fpace2/1000.;
5940 const G4int qrows = 50;
5941 const G4int qcols = 70;
5943 for(
G4int init_i = 0; init_i < qrows; init_i++) {
5944 for(
G4int init_j = 0; init_j < qcols; init_j++) {
5945 q[init_i][init_j] = 0.0;
5986 G4double aux=1.+(9.*xjj/4./qq/x13);
5988 G4double ee4=avol*xx+asur*(std::pow(xx,.666))+ac*x13+azer;
5989 G4double tota = ee1 + ee2 + ee3 + ee4;
5990 find = 939.55*xneu+938.77*zz - tota;
6002 const G4double fmp = 938.27231,fmn=939.56563,fml=1115.683;
6008 for(
G4int i=0;i<IMULTBU;i++){
6024 G4double gamma = std::sqrt(1.0 - v2 / (c*c));
6025 G4double avvmass = iz*fmp + (ia-iz-is)*fmn + is*fml +
eflmac(ia,iz,0,3);
6035 for(
G4int i=0;i<IEV_TAB;i++){
6051 G4double gamma = std::sqrt(1.0 - v2 / (c*c));
6052 G4double avvmass = iz*fmp + (ia-iz-is)*fmn + is*fml +
eflmac(ia,iz,0,3);
6065 G4double gamma = std::sqrt(1.0 - v2 / (c*c));
6139 return -1.0*std::abs(a);
6151 return -1*std::abs(a);
6160 fractpart = std::modf(number, &intpart);
6165 if(fractpart < 0.5) {
6166 return G4int(std::floor(number));
6169 return G4int(std::ceil(number));
6173 if(fractpart < -0.5) {
6174 return G4int(std::floor(number));
6177 return G4int(std::ceil(number));
6181 return G4int(std::floor(number));
6190 mylocaltime = localtime(&mytime);
6193 return(mylocaltime->tm_hour*60*60 + mylocaltime->tm_min*60 + mylocaltime->tm_sec);
6221 if(x-std::floor(x) <= std::ceil(x)-x)
6232 if(x-std::floor(x) <= std::ceil(x)-x)
6233 value =
G4int(std::floor(x));
6235 value =
G4int(std::ceil(x));
6242 if(x-std::floor(x) <= std::ceil(x)-x)
6243 return G4int(std::floor(x));
6245 return G4int(std::ceil(x));
6250 if(a < b && a < c) {
6253 if(b < a && b < c) {
6256 if(c < a && c < b) {
6276 G4int IZPART,IAPART,NMOTHER;
6279 G4double INT1,INT2,INT3,AKONST,EARG,R0,MPARTNER;
6282 G4double PAR_A1=0.,PAR_B1=0.,FACT=1.;
6297 NMOTHER =
idnint(AMOTHER-ZMOTHER);
6306 HBAR = 6.582122e-22;
6311 APARTNER = AMOTHER - APART;
6312 MPARTNER = APARTNER * 931.49 /
C2;
6315 if(IAPART==1&&IZPART==0){
6317 MPART = 939.56 /
C2;
6318 if(idlamb0==1)MPART = 1115.683 /
C2;
6320 if(IAPART==1&&IZPART==1){
6322 MPART = 938.27 /
C2;
6325 if(IAPART==2&&IZPART==0){
6327 MPART = 2.*939.56 /
C2;
6329 if(IAPART==2&&IZPART==1){
6331 MPART = 1876.10 /
C2;
6333 if(IAPART==3&&IZPART==1){
6335 MPART = 2809.39 /
C2;
6337 if(IAPART==3&&IZPART==2){
6339 MPART = 2809.37 /
C2;
6341 if(IAPART==4&&IZPART==2){
6343 MPART = 3728.35 /
C2;
6347 MPART = APART * 931.49 /
C2;
6357 MU = MPARTNER * MPART / (MPARTNER + MPART);
6361 RGEOM = R0 * (std::pow(APART,1.0/3.0)+std::pow(AMOTHER-APART,1.0/3.0));
6364 AKONST = HBAR*std::sqrt(1.0 / MU);
6367 BKONST = MPART / ( PI * PI * HBAR * HBAR);
6371 INT1 = 2.0 * std::pow(TEMP,3.) / (2.0 * TEMP +
B);
6373 ARG = std::sqrt(
B/TEMP);
6374 EARG = (
erf(ARG) - 1.0);
6375 if(std::abs(EARG)<1.e-9) EARG = 0.0;
6377 INT2 = 0.5 * std::sqrt(PI) * std::pow(TEMP,3.0/2.0);
6380 if(AEXP>700.0) AEXP = 700.0;
6381 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;
6382 if(INT2<0.0) INT2 = 0.0;
6385 if(EARG==0.0) INT2 = 0.0;
6388 INT3 = 2.0*TEMP*TEMP*TEMP / (2.0*TEMP*TEMP + 4.0*
B*TEMP +
B*
B);
6390 if(IZPART<-1.0&&ZMOTHER<151.0){
6394 fwidth = PI * BKONST * G * std::sqrt((RGEOM * RGEOM * INT1 + 2.0 * AKONST * RGEOM * INT2 + AKONST * AKONST * INT3) * RGEOM * RGEOM * INT1);
6397 fwidth = PI * BKONST * G *(RGEOM * RGEOM * INT1 + 2.0 * AKONST * RGEOM * INT2 + AKONST * AKONST * INT3);
6405 PAR_A1=std::exp(2.302585*0.2083*std::exp(-0.01548472*AMOTHER))-0.05;
6406 PAR_B1 = 0.59939389 + 0.00915657 * AMOTHER;
6408 if(AMOTHER>154.0&&AMOTHER<195.0){
6409 PAR_A1=1.0086961-8.629e-5*AMOTHER;
6410 PAR_B1 = 1.5329331 + 0.00302074 * AMOTHER;
6412 if(AMOTHER>194.0&&AMOTHER<208.0){
6413 PAR_A1=9.8356347-0.09294663*AMOTHER+2.441e-4*AMOTHER*AMOTHER;
6414 PAR_B1 = 7.7701987 - 0.02897401 * AMOTHER;
6416 if(AMOTHER>207.0&&AMOTHER<228.0){
6417 PAR_A1=15.107385-0.12414415*AMOTHER+2.7222e-4*AMOTHER*AMOTHER;
6418 PAR_B1=-64.078009+0.56813179*AMOTHER-0.00121078*AMOTHER*AMOTHER;
6421 if(
mod(NMOTHER,2)==0&&NMOTHER>147.){
6422 PAR_A1 = 2.0*(0.9389118 + 6.4559e-5 * AMOTHER);
6424 if(
mod(NMOTHER,2)==1)PAR_A1 = 3.0*(0.9389118 + 6.4559e-5 * AMOTHER);
6426 PAR_B1 = 2.1507177 + 0.00146119 * AMOTHER;
6432 FACT = std::exp((2.302585*PAR_A1*std::exp(-PAR_B1*(EXC-SB))));
6433 if(FACT<1.0) FACT = 1.0;
6434 if(IZPART<-1.&&ZMOTHER<151.0){
6436 fwidth = fwidth / std::sqrt(FACT);
6438 fwidth = fwidth / FACT;
6443 std::cout <<
"LOOK IN PARTICLE_WIDTH!" << std::endl;
6444 std::cout <<
"ACN,APART :"<< AMOTHER << APART << std::endl;
6445 std::cout <<
"EXC,TEMP,B,SB :" << EXC <<
" " << TEMP <<
" " <<
B <<
" " << SB << std::endl;
6446 std::cout <<
"INTi, i=1-3 :" << INT1 <<
" " << INT2 <<
" " << INT3 << std::endl;
6447 std::cout <<
" " << std::endl;
6468 HO = 197.3287 * omega;
6473 fpen=std::pow(10.0,4.e-4*std::pow(T/(HO*HO*std::pow(MU,0.25)),-4.3/2.3026));
6491 (*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;
6493 (*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;
6537 G4double DEFO_INIT,OMEGA,HOMEGA,OMEGA_GS,HOMEGA_GS,K1,MFCD;
6538 G4double BET1,XACT,SIGMA_SQR,W_EXP,XB,NORM,SIGMA_SQR_INF,W_INFIN,W;
6539 G4double FUNC_TRANS,LOG_SLOPE_INF,LOG_SLOPE_ABS;
6546 fomega_gs(AF,ZF,&K1,&OMEGA_GS,&HOMEGA_GS);
6550 if((bet*bet)>4.0*OMEGA_GS*OMEGA_GS){
6551 BET1=std::sqrt(bet*bet-4.0*OMEGA_GS*OMEGA_GS);
6556 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)));
6559 XACT = DEFO_INIT *std::exp(-0.5*(bet-BET1)*1.e21*(
TIME-T_0));
6564 BET1=std::sqrt(4.0*OMEGA_GS*OMEGA_GS-bet*bet);
6565 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));
6566 XACT = DEFO_INIT*std::cos(0.5*BET1*1.e21*(
TIME-T_0))*std::exp(-bet*1.e21*(
TIME-T_0));
6572 XB = 7./3.*
Y-938./765.*
Y*
Y+9.499768*
Y*
Y*
Y-8.050944*
Y*
Y*
Y*
Y;
6577 NORM = 1./std::sqrt(2.*PI*SIGMA_SQR);
6579 W_EXP = -1.*(XB - XACT)*(XB - XACT)/(2.0 * SIGMA_SQR);
6580 if(W_EXP<(-708.0) ) W_EXP = -708.0;
6581 W = NORM * std::exp( W_EXP ) * FT / (K1 * SIGMA_SQR);
6588 SIGMA_SQR_INF = FT/K1;
6589 W_EXP = -XB*XB/(2.0 * SIGMA_SQR_INF);
6590 if(W_EXP<(-708.0))W_EXP = -708.0;
6591 W_INFIN = std::exp(W_EXP)/std::sqrt(2.0*PI*SIGMA_SQR_INF);
6592 FUNC_TRANS = W / W_INFIN;
6597 LOG_SLOPE_INF =
cram(bet,HOMEGA)*bet*MFCD*OMEGA/FT;
6598 LOG_SLOPE_ABS = (XB-XACT)/SIGMA_SQR-XB/SIGMA_SQR_INF+
cram(bet,HOMEGA)*bet*MFCD*OMEGA/FT;
6600 FUNC_TRANS = FUNC_TRANS * LOG_SLOPE_ABS/LOG_SLOPE_INF;
6606void 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)
6665 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;
6683 if(BET*BET>4.0*OMEGA*OMEGA){
6690 t_0 = BET*1.e21*HBAR*HBAR/(4.*HOMEGA*FT)/16.;
6693 if(((2.*FT-HOMEGA/16.)>0.000001) && BET>0.0){
6698 t_0 = (std::log(2.*FT/(2.*FT-HOMEGA/16.)))/(BET*1.e21);
6710 STEP_LENGTH = 1.5*TAUF/50.;
6715 BEGIN_TIME = TSUM + t_0;
6717 if(BEGIN_TIME<0.0) std::cout <<
"CURRENT TIME < 0" << BEGIN_TIME << std::endl;
6719 if(BEGIN_TIME<1.50*TAUF){
6720 LOC_TIME_BEGIN = BEGIN_TIME;
6722 while((LOC_TIME_BEGIN<1.5*TAUF)&&fchoice==0){
6724 LOC_TIME_END = LOC_TIME_BEGIN + STEP_LENGTH;
6727 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;
6729 fGF_LOC = fGF_LOC * GF;
6739 LAMBDA = 1.0/TS1 + 1.0/TS2;
6745 REAC_PROB = std::exp(-1.0*STEP_LENGTH*LAMBDA);
6750 FISS_PROB = fGF_LOC / (fGF_LOC+GP);
6761 LOC_TIME_BEGIN = LOC_TIME_END;
6764 fT_LAPSE = LOC_TIME_END - BEGIN_TIME;
6771 FISS_PROB = GF / (GF+GP);
6781 LAMBDA = 1./TS1 + 1./TS2;
6803 (*T_LAPSE)=fT_LAPSE;
6816 G4double MFCD,OMEGA,HOMEGA1,HOMEGA2=0.,GFTUN;
6817 G4double E1,E2,EXP_FACT,CORR_FUNCT,FACT1,FACT2,FACT3;
6825 if(
mod(IN,2)==0&&
mod(IZ,2)==0){
6828 EE = EE - 12.0/std::sqrt(
A);
6832 if(
mod(IN,2)==1&&
mod(IZ,2)==1){
6836 if(
mod(IN,2)==1&&
mod(IZ,2)==0){
6840 if(
mod(IN,2)==0&&
mod(IZ,2)==1){
6844 E1 = EF + HOMEGA1/2.0/PI*std::log(HOMEGA1*(2.0*PI+HOMEGA2)/4.0/PI/PI);
6846 E2 = EF + HOMEGA2/(2.0*PI)*std::log(1.0+2.0*PI/HOMEGA2);
6853 EXP_FACT = (EE-EF)/(HOMEGA2/(2.0*PI));
6854 if(EXP_FACT>700.0) EXP_FACT = 700.0;
6855 CORR_FUNCT = HOMEGA1 * (1.0-1.0/(1.0+std::exp(EXP_FACT)));
6856 if(
mod(IN,2)==0&&
mod(IZ,2)==0){
6857 CORR_FUNCT = HOMEGA1 * (1.0-1.0/(1.0+std::exp(EXP_FACT)));
6860 FACT1 = HOMEGA1/(2.0*PI*TEMP+HOMEGA1);
6861 FACT2 = (2.0*PI/(2.0*PI+HOMEGA2)-HOMEGA1*(2.0*PI+HOMEGA2)/4.0/PI/PI)/(E2-E1);
6862 FACT3 = HOMEGA2/(2.0*PI*TEMP-HOMEGA2);
6865 GFTUN = FACT1*(std::exp(EE/TEMP)*std::exp(2.0*PI*(EE-EF)/HOMEGA1)-std::exp(-2.0*PI*EF/HOMEGA1));
6868 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));
6870 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));
6873 GFTUN = GFTUN/std::exp(EE/TEMP)*DENSF*ENH_FACT/DENSG/2.0/PI;
6874 GFTUN = GFTUN * CORR_FUNCT;
6879void 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)
6882 G4double FNORM,MASS_ASYM_SADD_B,FP_PER,FP_PAR,SIG_PER_SP,SIG_PAR_SP;
6883 G4double Z2OVERA,ftemp,fgf,DENSF,ECOR,EROT,qr;
6884 G4double DCR,UCR,ENH_FACTA,ENH_FACTB,ENH_FACT,PONFE;
6889 Z2OVERA = ZPRF * ZPRF /
A;
6892 if((ZPRF<=55.0) || (FF_ALLOWED==0)){
6902 densniv(
A,ZPRF,EE,EF,&DENSF,0.0,BS,BK,&ftemp,OPTSHP,0,
Y,&ECOR,JPR,1,&qr);
6905 fgf= DENSF/DENSG/PI/2.0*ftemp;
6917 FNORM = 1.2*1.2 * 931.49 * 1.e-2 / (9.0 * 6.582122*6.582122);
6920 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));
6925 if(Z2OVERA<=30.0) FP_PER = 6.50;
6928 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));
6929 if(FP_PAR<0.0) FP_PAR = 0.0;
6931 EROT = JPR * JPR / (2.0 * std::sqrt(FP_PAR*FP_PAR + FP_PER*FP_PER));
6932 if(IEROT==1) EROT = 0.0;
6935 SIG_PER_SP = std::sqrt(FP_PER * ftemp);
6937 if(SIG_PER_SP<1.0) SIG_PER_SP = 1.0;
6940 SIG_PAR_SP = std::sqrt(FP_PAR * ftemp);
6944 MASS_ASYM_SADD_B = 2.0;
6946 MASS_ASYM_SADD_B = 1.0;
6950 if(Z2OVERA>35.&&Z2OVERA<=(110.*110./298.0)){
6952 ENH_FACTA = std::sqrt(8.0*PI) * SIG_PER_SP*SIG_PER_SP * SIG_PAR_SP;
6954 ENH_FACTB = MASS_ASYM_SADD_B * SIG_PER_SP*SIG_PER_SP;
6956 ENH_FACT = ENH_FACTA * ENH_FACTB / (ENH_FACTA + ENH_FACTB);
6960 ENH_FACT = MASS_ASYM_SADD_B*SIG_PER_SP*SIG_PER_SP;
6963 ENH_FACT = std::sqrt(8.0*PI) * SIG_PER_SP*SIG_PER_SP* SIG_PAR_SP;
6968 PONFE = (ECOR-UCR-EROT)/DCR;
6969 if(PONFE>700.) PONFE = 700.0;
6971 ENH_FACT = 1.0/(1.0+std::exp(PONFE))*ENH_FACT+1.0;
6973 if(ENH_FACT<1.0)ENH_FACT = 1.0;
6974 fgf= DENSF/DENSG/PI/2.0*ftemp*ENH_FACT;
6978 fgf=
tunnelling(
A,ZPRF,
Y,EE,EF,ftemp,DENSG,DENSF,ENH_FACT);
6990 G4double AFRAGMENT,S4FINAL,ALEVDENS;
6991 G4double THETA_MOTHER,THETA_ORBITAL;
7006 if (EEFINAL<=0.01) EEFINAL = 0.01;
7007 AFRAGMENT = AMOTHER - ADAUGHTER;
7008 ALEVDENS = 0.073*AMOTHER + 0.095*std::pow(AMOTHER,2.0/3.0);
7009 S4FINAL = ALEVDENS * EEFINAL;
7010 if(S4FINAL <= 0.0 || S4FINAL > 100000.){
7011 std::cout<<
"S4FINAL:" << S4FINAL << ALEVDENS << EEFINAL <<
idnint(AMOTHER) <<
idnint(AFRAGMENT) << std::endl;
7013 THETA_MOTHER = 0.0111 * std::pow(AMOTHER,1.66667);
7014 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);
7016 *LORBITAL = -1.* THETA_ORBITAL * (LMOTHER / THETA_MOTHER + std::sqrt(S4FINAL) /(ALEVDENS*LMOTHER));
7018 *SIGMA_LORBITAL = std::sqrt(std::sqrt(S4FINAL) * THETA_ORBITAL / ALEVDENS);
7046 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;
7047 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;
7048 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;
7049 G4double PI = 3.141592653589793238;
7051 G4double SDWprevious=0,SUMDW_TOT=0,SUM_Z=0,X=0,SUMDW_N_TOT=0,XX=0;
7059 for (
G4int ia = 0; ia < 98; ia++)
7060 for (
G4int ib = 0; ib < 251; ib++) {
7061 BBIMF[ia][ib] = 0.0;
7062 SSBIMF[ia][ib] = 0.0;
7066 IZIMFMAX =
idnint(ZCN / 2.0);
7069 std::cout <<
"CHARGE_IMF line 46" << std::endl;
7070 std::cout <<
"Problem: IZIMFMAX < 3 " << std::endl;
7071 std::cout <<
"ZCN,IZIMFMAX," << ZCN <<
"," << IZIMFMAX << std::endl;
7079 bsbkbc(ACN,ZCN,&BS,&BK,&BC);
7081 densniv(ACN,ZCN,EE,0.0,&DENSCN,BSHELL,BS,BK,&TEMPCN,0,0,DEFBET,&ECOR,JPRF,0,&QRCN);
7083 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);
7084 EROTCN = JPRF * JPRF * 197.328 * 197.328 /(2. * IINERT);
7086 for(IZIMF=3;IZIMF<=IZIMFMAX;IZIMF++){
7095 INIMFMI =
max(1,INIMFMI-2);
7098 INCN =
idnint(ACN) - IZCN;
7102 INMI =
max(1,INMI-2);
7103 INMIN =
max(INIMFMI,INCN-INMA);
7104 INNMAX =
min(INIMFMA,INCN-INMI);
7106 ILIMMAX =
max(INNMAX,INMIN);
7109 for(
G4int INIMF=INMIN;INIMF<=ILIMMAX;INIMF++){
7110 IAIMF = IZIMF + INIMF;
7111 DW[IZIMF][IAIMF] = 0.0;
7112 AIMF_1 = 1.0*(IAIMF);
7115 mglms(ACN-AIMF_1,ZCN-ZIMF_1,OPTSHPIMF,&MARES);
7116 mglms(AIMF_1,ZIMF_1,OPTSHPIMF,&MAIMF);
7117 mglms(ACN,ZCN,OPTSHPIMF,&MAZ);
7121 SSBIMF[IZIMF][IAIMF] = 1.e37;
7124 SSBIMF[IZIMF][IAIMF] = MAIMF + MARES - MAZ + fBIMF;
7125 BBIMF[IZIMF][IAIMF] = fBIMF;
7131 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));
7133 EROT = JPRF * JPRF * 197.328 * 197.328 /(2.0 * IINERT);
7136 if (EE<(SSBIMF[IZIMF][IAIMF]+EROT) || DENSCN<=0.0){
7145 densniv(ACN,ZCN,EE,SSBIMF[IZIMF][IAIMF],&DENSIMF,0.0,BSIMF,1.0,&fTIMF,0,0,DEFBETIMF,&ECOR,JPRF,2,&QR);
7146 IMFARG = (SSBIMF[IZIMF][IAIMF]+EROTCN-EROT)/fTIMF;
7147 if(IMFARG>200.0) IMFARG = 200.0;
7149 WIDTH1 =
width(ACN,ZCN,AIMF_1,ZIMF_1,fTIMF,fBIMF,SSBIMF[IZIMF][IAIMF],EE-EROT);
7151 WIDTH_IMF = WIDTH1 * std::exp(-IMFARG) * QR / QRCN;
7154 std::cout <<
"GAMMA_IMF=0 -> LOOK IN GAMMA_IMF CALCULATIONS!" << std::endl;
7155 std::cout <<
"ACN,ZCN,AIMF,ZIMF:" <<
idnint(ACN) <<
"," <<
idnint(ZCN) <<
"," <<
idnint(AIMF_1) <<
"," <<
idnint(ZIMF_1) << std::endl;
7156 std::cout <<
"SSBIMF,TIMF :" << SSBIMF[IZIMF][IAIMF] <<
"," << fTIMF << std::endl;
7157 std::cout <<
"DEXP(-IMFARG) = " << std::exp(-IMFARG) << std::endl;
7158 std::cout <<
"WIDTH1 =" << WIDTH1 << std::endl;
7162 SDW[IZIMF] = SDW[IZIMF] + WIDTH_IMF;
7164 DW[IZIMF][IAIMF] = WIDTH_IMF;
7172 SDWprevious = 1.e20;
7175 for(
G4int III_ZIMF=3;III_ZIMF<=IZIMFMAX;III_ZIMF++){
7177 if(SDW[III_ZIMF]==0.0){
7178 IZSTOP = III_ZIMF - 1;
7182 if(SDW[III_ZIMF]>SDWprevious){
7183 IZSTOP = III_ZIMF - 1;
7186 SDWprevious = SDW[III_ZIMF];
7198 A1PAR = std::log10(SDW[IZSTOP]/SDW[IZSTOP-2])/std::log10((1.0*IZSTOP)/(1.0*IZSTOP-2.0));
7199 A2PAR = std::log10(SDW[IZSTOP]) - A1PAR * std::log10(1.0*(IZSTOP));
7200 if(A2PAR>0.)A2PAR=-1.*A2PAR;
7201 if(A1PAR>0.)A1PAR=-1.*A1PAR;
7205 for(
G4int II_ZIMF = IZSTOP;II_ZIMF<=IZIMFMAX;II_ZIMF++){
7206 SDW[II_ZIMF] = std::pow(10.0,A2PAR) * std::pow(1.0*II_ZIMF,A1PAR);
7207 if(SDW[II_ZIMF]<0.0) SDW[II_ZIMF] = 0.0;
7214 for(
G4int I_ZIMF = 3;I_ZIMF<=IZIMFMAX;I_ZIMF++){
7215 SUMDW_TOT = SUMDW_TOT + SDW[I_ZIMF];
7218 std::cout <<
"*********************" << std::endl;
7219 std::cout <<
"IMF function" << std::endl;
7220 std::cout <<
"SUM of decay widths = " << SUMDW_TOT <<
" IZIMFMAX = " << IZIMFMAX << std::endl;
7221 std::cout <<
"IZSTOP = " << IZSTOP << std::endl;
7229 X =
haz(1)*SUMDW_TOT;
7236 for(
G4int IZ = 3;IZ<=IZIMFMAX;IZ++){
7237 SUM_Z = SUM_Z + SDW[IZ];
7250 INMINMEM =
max(1,INMINMEM-2);
7253 INMI =
max(1,INMI-2);
7256 INMINMEM =
max(INMINMEM,INCN-INMA);
7257 INMAXMEM =
min(INMAXMEM,INCN-INMI);
7259 INMAXMEM =
max(INMINMEM,INMAXMEM);
7263 for(
G4int IIINIMF = INMINMEM;IIINIMF<=INMAXMEM;IIINIMF++){
7264 IA = IZMEM + IIINIMF;
7265 if(IZMEM>=3&&IZMEM<=95&&IA>=4&&IA<=250){
7266 SUMDW_N_TOT = SUMDW_N_TOT + DW[IZMEM][IA];
7268 std::cout <<
"CHARGE IMF OUT OF RANGE" << IZMEM <<
", " << IA <<
", " <<
idnint(ACN) <<
", " <<
idnint(ZCN) <<
", " << TEMP << std::endl;
7272 XX =
haz(1)*SUMDW_N_TOT;
7275 for(
G4int IINIMF = INMINMEM;IINIMF<=INMAXMEM; IINIMF++){
7276 IIA = IZMEM + IINIMF;
7278 SUM_A = SUM_A + DW[IZMEM][IIA];
7287 NIMF = fAIMF - fZIMF;
7289 if((ACN-ZCN-NIMF)<=0.0 || (ZCN-fZIMF) <= 0.0){
7290 std::cout <<
"IMF Partner unstable:" << std::endl;
7291 std::cout <<
"System: Acn,Zcn,NCN:" << std::endl;
7292 std::cout <<
idnint(ACN) <<
", " <<
idnint(ZCN) <<
", " <<
idnint(ACN-ZCN) << std::endl;
7293 std::cout <<
"IMF: A,Z,N:" << std::endl;
7294 std::cout <<
idnint(fAIMF) <<
", " <<
idnint(fZIMF) <<
", " <<
idnint(fAIMF-fZIMF) << std::endl;
7295 std::cout <<
"Partner: A,Z,N:" << std::endl;
7296 std::cout <<
idnint(ACN-fAIMF) <<
", " <<
idnint(ZCN-fZIMF) <<
", " <<
idnint(ACN-ZCN-NIMF) << std::endl;
7297 std::cout <<
"----nmin,nmax" << INMINMEM <<
", " << INMAXMEM << std::endl;
7298 std::cout <<
"----- warning: Zimf=" << fZIMF <<
" Aimf=" << fAIMF << std::endl;
7299 std::cout <<
"----- look in subroutine IMF" << std::endl;
7300 std::cout <<
"ACN,ZCN,ZIMF,AIMF,temp,EE,JPRF::" << ACN <<
", " << ZCN <<
", " << fZIMF <<
", " << fAIMF <<
", " << TEMP <<
", " << EE <<
", " << JPRF << std::endl;
7301std::cout <<
"-IZSTOP,IZIMFMAX:" << IZSTOP <<
", " << IZIMFMAX << std::endl;
7302std::cout <<
"----X,SUM_Z,SUMDW_TOT:" << X <<
", " << SUM_Z <<
", " << SUMDW_TOT << std::endl;
7308 if(fZIMF>=ZCN || fAIMF>=ACN || fZIMF<=2 || fAIMF<=3){
7309 std::cout <<
"----nmin,nmax" << INMINMEM <<
", " << INMAXMEM << std::endl;
7310 std::cout <<
"----- warning: Zimf=" << fZIMF <<
" Aimf=" << fAIMF << std::endl;
7311 std::cout <<
"----- look in subroutine IMF" << std::endl;
7312 std::cout <<
"ACN,ZCN,ZIMF,AIMF,temp,EE,JPRF:" << ACN <<
", " << ZCN <<
", " << fZIMF <<
", " << fAIMF <<
", " << TEMP <<
", " << EE <<
", " << JPRF << std::endl;
7313std::cout <<
"-IZSTOP,IZIMFMAX:" << IZSTOP <<
", " << IZIMFMAX << std::endl;
7314std::cout <<
"----X,SUM_Z,SUMDW_TOT:" << X <<
", " << SUM_Z <<
", " << SUMDW_TOT << std::endl;
7315for(
int III_ZIMF=3;III_ZIMF<=IZIMFMAX;III_ZIMF++)
7316 std::cout <<
"-**Z,SDW:" << III_ZIMF <<
", " << SDW[III_ZIMF] << std::endl;
7326 if((ZCN-fZIMF)<=0.0)std::cout <<
"CHARGE_IMF ZIMF > ZCN" << std::endl;
7327 if((ACN-fAIMF)<=0.0)std::cout <<
"CHARGE_IMF AIMF > ACN" << std::endl;
7332 EEDAUG = (EE - fSBIMF) * (ACN - fAIMF) / ACN;
7333 bsbkbc(ACN - fAIMF,ZCN-fZIMF,&BS,&BK,&BC);
7334 densniv(ACN-fAIMF,ZCN-fZIMF,EEDAUG,0.0,&DENSIMF,BSHELL,BS,BK,&fTIMF,0,0,DEFBET,&ECOR,0.0,0,&QR);
7337 std::cout <<
"----- warning: EE=" << EE <<
"," <<
" S+Bimf=" << fSBIMF << std::endl;
7338 std::cout <<
"----- look in subroutine IMF" << std::endl;
7339 std::cout <<
"IMF will be resampled" << std::endl;
7353G4int VISOSTAB[191][2]={
7464 *nmin = VISOSTAB[z-1][0];
7465 *nmax = VISOSTAB[z-1][1];
7481 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;
7482 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;
7483 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;
7485 G4double xcv=0.,ycv=0.,zcv=0.,VXOUT=0.,VYOUT=0.,VZOUT=0.;
7487 G4double jprfn=0.0, jprfp=0.0, jprfd=0.0, jprft=0.0, jprfhe=0.0, jprfa=0.0, jprflamb0=0.0;
7488 G4double ctet1 = 0.0, stet1 = 0.0, phi1 = 0.0;
7491 G4int itest = 0, sortie=0;
7497 G4double time,tauf,tau0,
a0,a1,emin,ts1,tsum=0.;
7498 G4int inttype=0,inum=0,gammadecay = 0, flamb0decay = 0;
7504 G4int NbLam0= (*NbLam0_par);
7508 const G4double mu2 = 931.494*931.494;
7528 a0 = 0.66482503 - 3.4678935 * std::exp(-0.0104002*ee);
7529 a1 = 5.6846e-04 + 0.00574515 * std::exp(-0.01114307*ee);
7530 tauf = (
a0 + a1 * zf*zf/std::pow(af,0.3333333)) * tau0;
7533 direct(zf,af,ee,0.,&probp,&probd,&probt,&probn,&probhe,&proba,&probg,&probimf,&probf,&problamb0,&ptotl,
7534 &sn,&sbp,&sbd,&sbt,&sbhe,&sba,&slamb0,
7535 &ecn,&ecp,&ecd,&ect,&eche,&eca,&ecg,&eclamb0,
7536 &
bp,&bd,&bt,&bhe,&ba,&
sp,&sd,&st,&she,&sa,&ef,&ts1,inttype,inum,itest,&sortie,&tcn,
7537 &jprfn, &jprfp, &jprfd, &jprft, &jprfhe, &jprfa, &jprflamb0, &tsum, NbLam0);
7541 if(ptotl<=0.)
goto post100;
7545 if(emin>1e30)std::cout <<
"ERROR AT THE EXIT OF EVAPORA,E>1.D30,AF" << std::endl;
7552 pc = std::sqrt(std::pow((1.0 + ecn/9.3956e2),2.) - 1.0) * 9.3956e2;
7556 else if(probp!=0.0){
7560 pc = std::sqrt(std::pow((1.0 + ecp/9.3827e2),2.) - 1.0) * 9.3827e2;
7564 else if(probd!=0.0){
7568 pc = std::sqrt(std::pow((1.0 + ecd/1.875358e3),2) - 1.0) * 1.875358e3;
7572 else if(probt!=0.0){
7576 pc = std::sqrt(std::pow((1.0 + ect/2.80828e3),2) - 1.0) * 2.80828e3;
7580 else if(probhe!=0.0){
7583 epsiln = she + eche;
7584 pc = std::sqrt(std::pow((1.0 + eche/2.80826e3),2) - 1.0) * 2.80826e3;
7588 else{
if(proba!=0.0){
7592 pc = std::sqrt(std::pow((1.0 + eca/3.72834e3),2) - 1.0) * 3.72834e3;
7615 pc = std::sqrt(std::pow((1.0 + eca/3.72834e3),2) - 1.0) * 3.72834e3;
7619 else if (x < proba+probhe) {
7623 epsiln = she + eche;
7624 pc = std::sqrt(std::pow((1.0 + eche/2.80826e3),2) - 1.0) * 2.80826e3;
7628 else if (x < proba+probhe+probt) {
7633 pc = std::sqrt(std::pow((1.0 + ect/2.80828e3),2) - 1.0) * 2.80828e3;
7637 else if (x < proba+probhe+probt+probd) {
7642 pc = std::sqrt(std::pow((1.0 + ecd/1.875358e3),2) - 1.0) * 1.875358e3;
7646 else if (x < proba+probhe+probt+probd+probp) {
7651 pc = std::sqrt(std::pow((1.0 + ecp/9.3827e2),2) - 1.0) * 9.3827e2;
7655 else if (x < proba+probhe+probt+probd+probp+probn) {
7660 pc = std::sqrt(std::pow((1.0 + ecn/9.3956e2),2.) - 1.0) * 9.3956e2;
7664 else if (x < proba+probhe+probt+probd+probp+probn+problamb0) {
7668 epsiln = slamb0 + eclamb0;
7669 pc = std::sqrt(std::pow((1.0 + (eclamb0)/11.1568e2),2.) - 1.0) * 11.1568e2;
7675 else if (x < proba+probhe+probt+probd+probp+probn+problamb0+probg) {
7683 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){
7694 if(gammadecay==1 && ee<=0.01+epsiln){
7703 if(ee<=0.01) ee = 0.010;
7705 if(af<2.5)
goto post100;
7721 ctet1 = 2.0*rnd - 1.0;
7722 stet1 = std::sqrt(1.0 - std::pow(ctet1,2));
7724 phi1 = rnd*2.0*3.141592654;
7725 xcv = stet1*std::cos(phi1);
7726 ycv = stet1*std::sin(phi1);
7731 G4double ETOT_LP = std::sqrt(
pc*
pc + amoins*amoins * mu2);
7732 if(flamb0decay==1)ETOT_LP = std::sqrt(
pc*
pc + 1115.683*1115.683);
7745 &VXOUT,&VYOUT,&VZOUT);
7755 G4double gamma = 1.0/std::sqrt(1.0 - v2 / (c*c));
7756 G4double etot_lp = amoins*mu * gamma;
7766 pteva = std::sqrt(pxeva*pxeva + pyeva*pyeva);
7768 etot = std::sqrt ( pleva*pleva + pteva*pteva + af*af * mu2 );
7769 vx_eva = c * pxeva / etot;
7770 vy_eva = c * pyeva / etot;
7771 vz_eva = c * pleva / etot;
7775 if(time<tauf)
goto post10;
7781 *E_scission_post = ee;
7782 *NbLam0_par = NbLam0;
7787 if(
A<1.)
return (1.*H)/
A*(10.68*
A-21.27*std::pow(
A,2./3.))*10.;
7788 return (1.*H)/
A*(10.68*
A-21.27*std::pow(
A,2./3.));
7793 if(
A<1.)
return 1.e38;
7797 if (
Z == 1 &&
A == 4)
7799 else if (
Z == 2 &&
A == 4)
7801 else if (
Z == 2 &&
A == 5)
7803 else if (
Z == 2 &&
A == 6)
7805 else if (
Z == 2 &&
A == 7)
7807 else if (
Z == 2 &&
A == 8)
7809 else if (
Z == 3 &&
A == 6)
7811 else if (
Z == 3 &&
A == 7)
7813 else if (
Z == 3 &&
A == 8)
7815 else if (
Z == 3 &&
A == 9)
7817 else if (
Z == 4 &&
A == 7)
7819 else if (
Z == 4 &&
A == 8)
7821 else if (
Z == 4 &&
A == 9)
7823 else if (
Z == 4 &&
A == 10)
7825 else if (
Z == 5 &&
A == 9)
7827 else if (
Z == 5 &&
A == 10)
7829 else if (
Z == 5 &&
A == 11)
7831 else if (
Z == 5 &&
A == 12)
7833 else if (
Z == 6 &&
A == 12)
7835 else if (
Z == 6 &&
A == 13)
7837 else if (
Z == 6 &&
A == 14)
7839 else if (
Z == 7 &&
A == 14)
7841 else if (
Z == 7 &&
A == 15)
7843 else if (
Z == 8 &&
A == 16)
7845 else if (
Z == 8 &&
A == 17)
7847 else if (
Z == 14 &&
A == 28)
7849 else if (
Z == 39 &&
A == 89)
7851 else if (
Z == 57 &&
A == 139)
7853 else if (
Z == 82 &&
A == 208)
7865 if(
A<2 ||
Z<2)
return 0.;
7875 if(
mod(N,2) == 1 &&
mod(
Z,2) == 1)
D = -12./std::sqrt(
A);
7876 if(
mod(N,2) == 0 &&
mod(
Z,2) == 0)
D = 12./std::sqrt(
A);
7878 G4double deltanew = (1.-std::exp(-1.*
A/c))*
D;
7880 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));
7884void 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)
7893 e =
dmin1(SBHE,SN,e);
7894 e =
dmin1(SBHE,SBA,e);
7915 *ECP = (-1.0)*SP +
BP;
7932 *ECD = (-1.0)*SD + BD;
7949 *ECT = (-1.0)*ST + BT;
7966 *ECHE= (-1.0)*SHE + BHE;
7984 *ECA = (-1.0)*SA + BA;
8066 Delta_U1_shell_max = -2.45;
8072 Delta_U2_shell = -2.45;
8082 G4double Fwidth_asymm1,Fwidth_asymm2,Fwidth_symm;
8084 Fwidth_asymm1 = 0.65;
8085 Fwidth_asymm2 = 0.65;
8127 Friction_factor = 1.0;
8130 G4double cN_asymm1_shell, cN_asymm2_shell;
8131 G4double gamma,gamma_heavy1,gamma_heavy2;
8147 G4double Sasymm1=0.,Sasymm2=0.,Ssymm=0.,Ysum=0.,Yasymm=0.;
8149 G4double wNasymm1_saddle, wNasymm2_saddle, wNsymm_saddle;
8150 G4double wNasymm2_scission, wNsymm_scission;
8151 G4double wNasymm1, wNasymm2, wNsymm;
8174 G4double A_levdens_heavy1,A_levdens_heavy2;
8178 G4double epsilon_1_saddle,epsilon0_1_saddle;
8179 G4double epsilon_2_saddle,epsilon0_2_saddle,epsilon_symm_saddle;
8184 G4double E_eff1_saddle,E_eff2_saddle;
8185 G4double Epot0_mode1_saddle,Epot0_mode2_saddle,Epot0_symm_saddle;
8186 G4double Epot_mode1_saddle,Epot_mode2_saddle,Epot_symm_saddle;
8187 G4double E_defo,E_defo1,E_defo2,E_scission_pre=0.,E_scission_post;
8193 G4double MassCurv_scis, MassCurv_sadd;
8195 G4double Nheavy1_shell,Nheavy2_shell;
8200 G4double Z_scission,N_scission,A_scission;
8202 G4double beta1gs=0.,beta2gs=0.,betags=0.;
8204 G4double DSN132,Delta_U1_shell,E_eff0_saddle;
8205 G4int NbLam0= (*NbLam0_par);
8210 cN_asymm1_shell = 0.700 * N/
Z;
8211 cN_asymm2_shell = 0.040 * N/
Z;
8215 DSN132 = Nheavy1_in - N/
Z * Zheavy1_in;
8216 Aheavy1 = Nheavy1_in + Zheavy1_in + 0.340 * DSN132;
8225 Delta_U1_shell = Delta_U1_shell_max + U1NZ_SLOPE * std::abs(DSN132);
8226 Delta_U1_shell =
min(0.,Delta_U1_shell);
8230 Nheavy1 = N/
A * Aheavy1;
8231 Aheavy2 = Nheavy2 *
A/N;
8235 A_levdens =
A / xLevdens;
8236 gamma = A_levdens / (0.40 * std::pow(
A,1.3333)) * FGAMMA;
8237 A_levdens_heavy1 = Aheavy1 / xLevdens;
8238 gamma_heavy1 = A_levdens_heavy1 / (0.40 * std::pow(Aheavy1,1.3333)) * FGAMMA * FGAMMA1;
8239 A_levdens_heavy2 = Aheavy2 / xLevdens;
8240 gamma_heavy2 = A_levdens_heavy2 / (0.40 * std::pow(Aheavy2,1.3333)) * FGAMMA;
8244 E_saddle_scission = (-24. + 0.02227 *
Z*
Z/std::pow(
A,0.33333))*Friction_factor;
8245 E_saddle_scission =
max( 0.0, E_saddle_scission );
8251 Z2_over_A_eff =
Z*
Z/
A;
8253 if( Z2_over_A_eff< 34.0 )
8254 MassCurv_scis = std::pow(10., -1.093364 + 0.082933 * Z2_over_A_eff - 0.0002602 * Z2_over_A_eff*Z2_over_A_eff);
8256 MassCurv_scis = std::pow(10., 3.053536 - 0.056477 * Z2_over_A_eff+ 0.0002454 * Z2_over_A_eff*Z2_over_A_eff );
8261 MassCurv_sadd = X_s2s * MassCurv_scis;
8263 cN_symm = 8.0 / std::pow(N,2.) * MassCurv_scis;
8264 cN_symm_sadd = 8.0 / std::pow(N,2.) * MassCurv_sadd;
8266 Nheavy1_shell = Nheavy1;
8269 Nheavy1_eff = (cN_symm_sadd*Nsymm + cN_asymm1_shell *
8270 Uwash(E/
A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1) *
8274 Uwash(E/
A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1));
8276 Nheavy1_eff = (cN_symm_sadd*Nsymm +
8277 cN_asymm1_shell*Nheavy1_shell)
8282 Nheavy2_NZ = Nheavy2;
8283 Nheavy2_shell = Nheavy2_NZ;
8285 Nheavy2_eff = (cN_symm_sadd*Nsymm +
8287 Uwash(E/
A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2) *
8291 Uwash(E/
A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2));
8293 Nheavy2_eff = (cN_symm_sadd*Nsymm +
8294 cN_asymm2_shell*Nheavy2_shell)
8298 Delta_U1 = Delta_U1_shell + (Nheavy1_shell - Nheavy1_eff)*(Nheavy1_shell - Nheavy1_eff) * cN_asymm1_shell;
8299 Delta_U1 =
min(Delta_U1,0.0);
8300 Delta_U2 = Delta_U2_shell + (Nheavy2_shell - Nheavy2_eff)*(Nheavy2_shell - Nheavy2_eff) * cN_asymm2_shell;
8301 Delta_U2 =
min(Delta_U2,0.0);
8305 Epot0_mode1_saddle = (Nheavy1_eff-Nsymm)*(Nheavy1_eff-Nsymm) * cN_symm_sadd;
8306 Epot0_mode2_saddle = (Nheavy2_eff-Nsymm)*(Nheavy2_eff-Nsymm) * cN_symm_sadd;
8307 Epot0_symm_saddle = 0.0;
8311 Epot_mode1_saddle = Epot0_mode1_saddle + Delta_U1;
8312 Epot_mode2_saddle = Epot0_mode2_saddle + Delta_U2;
8313 Epot_symm_saddle = Epot0_symm_saddle;
8316 dUeff =
min( Epot_mode1_saddle, Epot_mode2_saddle);
8317 dUeff =
min( dUeff, Epot_symm_saddle);
8318 dUeff = dUeff - Epot_symm_saddle;
8327 epsilon0_1_saddle = Eld - Epot0_mode1_saddle;
8328 epsilon0_2_saddle = Eld - Epot0_mode2_saddle;
8331 epsilon_1_saddle = Eld - Epot_mode1_saddle;
8332 epsilon_2_saddle = Eld - Epot_mode2_saddle;
8334 epsilon_symm_saddle = Eld - Epot_symm_saddle;
8337 eexc1_saddle = epsilon_1_saddle;
8338 eexc2_saddle = epsilon_2_saddle;
8341 EEXC_MAX =
max( eexc1_saddle, eexc2_saddle);
8342 EEXC_MAX =
max( EEXC_MAX, Eld);
8345 epsilon_1_scission = Eld + E_saddle_scission - Epot_mode1_saddle;
8346 epsilon_2_scission = Eld + E_saddle_scission - Epot_mode2_saddle;
8349 epsilon_symm_scission = Eld + E_saddle_scission - Epot_symm_saddle;
8352 E_eff1_saddle = epsilon0_1_saddle - Delta_U1 *
8353 Uwash(epsilon_1_saddle/
A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1);
8355 if( E_eff1_saddle < A_levdens * hbom1*hbom1)
8356 E_eff1_saddle = A_levdens * hbom1*hbom1;
8359 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_eff1_saddle) /
8361 Uwash(epsilon_1_saddle/
A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)+
8364 E_eff2_saddle = epsilon0_2_saddle -
8366 Uwash(epsilon_2_saddle/
A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2);
8368 if(E_eff2_saddle < A_levdens * hbom2*hbom2)
8369 E_eff2_saddle = A_levdens * hbom2*hbom2;
8372 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_eff2_saddle) /
8374 Uwash(epsilon_2_saddle/
A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)+
8377 E_eff0_saddle = epsilon_symm_saddle;
8378 if(E_eff0_saddle < A_levdens * hbom3*hbom3)
8379 E_eff0_saddle = A_levdens * hbom3*hbom3;
8382 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_eff0_saddle) /
8385 if(epsilon_symm_scission > 0.0 ){
8386 E_HELP =
max(E_saddle_scission,epsilon_symm_scission);
8388 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*(E_HELP)) /
8392 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_saddle_scission) /
8399 if( E_saddle_scission == 0.0 ){
8400 wNasymm1_scission = wNasymm1_saddle;
8401 wNasymm2_scission = wNasymm2_saddle;
8403 if( Nheavy1_eff > 75.0 ){
8404 wNasymm1_scission = std::sqrt(21.0)*N/
A;
8405 wNasymm2_scission =
max( 12.8 - 1.0 *(92.0 - Nheavy2_eff),1.0)*N/
A;
8408 wNasymm1_scission = wNasymm1_saddle;
8409 wNasymm2_scission = wNasymm2_saddle;
8413 wNasymm1_scission =
max( wNasymm1_scission, wNasymm1_saddle );
8414 wNasymm2_scission =
max( wNasymm2_scission, wNasymm2_saddle );
8416 wNasymm1 = wNasymm1_scission * Fwidth_asymm1;
8417 wNasymm2 = wNasymm2_scission * Fwidth_asymm2;
8418 wNsymm = wNsymm_scission * Fwidth_symm;
8421 Aheavy1_eff = Nheavy1_eff *
A/N;
8422 Aheavy2_eff = Nheavy2_eff *
A/N;
8424 A_levdens_heavy1 = Aheavy1_eff / xLevdens;
8425 A_levdens_heavy2 = Aheavy2_eff / xLevdens;
8426 gamma_heavy1 = A_levdens_heavy1 / (0.40 * std::pow(Aheavy1_eff,1.3333)) * FGAMMA * FGAMMA1;
8427 gamma_heavy2 = A_levdens_heavy2 / (0.40 * std::pow(Aheavy2_eff,1.3333)) * FGAMMA;
8429 if( epsilon_symm_saddle < A_levdens * hbom3*hbom3)
8430 Ssymm = 2.0 * std::sqrt(A_levdens*A_levdens * hbom3*hbom3) +
8431 (epsilon_symm_saddle - A_levdens * hbom3*hbom3)/hbom3;
8433 Ssymm = 2.0 * std::sqrt(A_levdens*epsilon_symm_saddle);
8437 if( epsilon0_1_saddle < A_levdens * hbom1*hbom1 )
8438 Ssymm_mode1 = 2.0 * std::sqrt(A_levdens*A_levdens * hbom1*hbom1) +
8439 (epsilon0_1_saddle - A_levdens * hbom1*hbom1)/hbom1;
8441 Ssymm_mode1 = 2.0 * std::sqrt( A_levdens*epsilon0_1_saddle );
8443 if( epsilon0_2_saddle < A_levdens * hbom2*hbom2 )
8444 Ssymm_mode2 = 2.0 * std::sqrt(A_levdens*A_levdens * hbom2*hbom2) +
8445 (epsilon0_2_saddle - A_levdens * hbom2*hbom2)/hbom2;
8447 Ssymm_mode2 = 2.0 * std::sqrt(A_levdens*epsilon0_2_saddle);
8450 if( epsilon0_1_saddle -
8452 Uwash(epsilon_1_saddle/
A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)
8453 < A_levdens * hbom1*hbom1 )
8454 Sasymm1 = 2.0 * std::sqrt( A_levdens*A_levdens * hbom1*hbom1 ) +
8455 (epsilon0_1_saddle - Delta_U1 *
8456 Uwash(epsilon_1_saddle/
A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)
8457 - A_levdens * hbom1*hbom1)/hbom1;
8459 Sasymm1 = 2.0 *std::sqrt( A_levdens*(epsilon0_1_saddle - Delta_U1 *
8460 Uwash(epsilon_1_saddle/
A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)));
8462 if( epsilon0_2_saddle -
8464 Uwash(epsilon_2_saddle/
A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)
8465 < A_levdens * hbom2*hbom2 )
8466 Sasymm2 = 2.0 * std::sqrt( A_levdens*A_levdens * hbom2*hbom2 ) +
8467 (epsilon0_1_saddle-Delta_U1 *
8468 Uwash(epsilon_2_saddle/
A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)
8469 - A_levdens * hbom2*hbom2)/hbom2;
8472 std::sqrt( A_levdens*(epsilon0_2_saddle - Delta_U2 *
8473 Uwash(epsilon_2_saddle/
A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)));
8475 Yasymm1 = ( std::exp(Sasymm1 - Ssymm) - std::exp(Ssymm_mode1 - Ssymm) ) *
8476 wNasymm1_saddle / wNsymm_saddle * 2.0;
8478 Yasymm2 = ( std::exp(Sasymm2 - Ssymm) - std::exp(Ssymm_mode2 - Ssymm) ) *
8479 wNasymm2_saddle / wNsymm_saddle * 2.0;
8481 Ysum = Ysymm + Yasymm1 + Yasymm2;
8484 Ysymm = Ysymm / Ysum;
8485 Yasymm1 = Yasymm1 / Ysum;
8486 Yasymm2 = Yasymm2 / Ysum;
8487 Yasymm = Yasymm1 + Yasymm2;
8493 if( (epsilon_symm_saddle < epsilon_1_saddle) &&
8494 (epsilon_symm_saddle < epsilon_2_saddle) )
8497 if( epsilon_1_saddle < epsilon_2_saddle )
8505 r_e_o = std::pow(10.0,-0.0170 * (E_saddle_scission + Eld)*(E_saddle_scission + Eld));
8523 if( rmode < Yasymm1 )
8526 if( (rmode > Yasymm1) && (rmode < Yasymm) )
8535 N1mean = Nheavy1_eff;
8539 N1mean = Nheavy2_eff;
8555 while( N1r < 5.0 || N2r < 5.0 ){
8574 E_scission_pre =
max( epsilon_1_scission, 1.0 );
8576 if( N1mean > N*0.50 ){
8586 E_scission_pre =
max( epsilon_2_scission, 1.0 );
8587 if( N1mean > N*0.50 ){
8588 beta1 = (N1r - 92.0) * 0.030 + 0.60;
8593 beta1 =
max(beta1,beta1gs);
8594 beta2 = 1.0 - beta1;
8595 beta2 =
max(beta2,beta2gs);
8601 beta2 = (N2r -92.0) * 0.030 + 0.60;
8602 beta2 =
max(beta2,beta2gs);
8603 beta1 = 1.0 - beta2;
8604 beta1 =
max(beta1,beta1gs);
8619 beta =
max(0.177963+0.0153241*Zsymm-1.62037e-4*Zsymm*Zsymm,betags);
8620 beta1 =
max(0.177963+0.0153241*Z1UCD-1.62037e-4*Z1UCD*Z1UCD,beta1gs);
8621 beta2 =
max(0.177963+0.0153241*Z2UCD-1.62037e-4*Z2UCD*Z2UCD,beta2gs);
8623 E_asym =
frldm( Z1UCD, N1r, beta1 ) +
8624 frldm( Z2UCD, N2r, beta2 ) +
8625 ecoul( Z1UCD, N1r, beta1, Z2UCD, N2r, beta2, 2.0 ) -
8628 E_scission_pre =
max( epsilon_symm_scission - E_asym, 1. );
8637 if(E_scission_pre>5. && NbLam0<1){
8639 &A_scission,&Z_scission,vx_eva_sc,vy_eva_sc,vz_eva_sc,&NbLam0);
8640 N_scission = A_scission - Z_scission;
8644 E_scission_post = E_scission_pre;
8645 N_scission = A_scission - Z_scission;
8651 N1r = N1r * N_scission / N;
8652 N2r = N2r * N_scission / N;
8653 Z1UCD = Z1UCD * Z_scission /
Z;
8654 Z2UCD = Z2UCD * Z_scission /
Z;
8692 CZ = (
frldm( Z1UCD-1.0, N1r+1.0, beta1 ) +
8693 frldm( Z2UCD+1.0, N2r-1.0, beta2 ) +
8694 frldm( Z1UCD+1.0, N1r-1.0, beta1 ) +
8695 frldm( Z2UCD-1.0, N2r+1.0, beta2 ) +
8696 ecoul( Z1UCD-1.0, N1r+1.0, beta1,
8697 Z2UCD+1.0, N2r-1.0, beta2, 2.0) +
8698 ecoul( Z1UCD+1.0, N1r-1.0, beta1,
8699 Z2UCD-1.0, N2r+1.0, beta2, 2.0) -
8700 2.0*
ecoul( Z1UCD, N1r, beta1, Z2UCD, N2r, beta2, 2.0) -
8701 2.0*
frldm( Z1UCD, N1r, beta1 ) -
8702 2.0*
frldm( Z2UCD, N2r, beta2) ) * 0.50;
8704 if(1.0/A_levdens*E_scission_post < 0.0)
8705 std::cout <<
"DSQRT 1 < 0" << A_levdens <<
" " << E_scission_post << std::endl;
8707 if(0.50 * std::sqrt(1.0/A_levdens*E_scission_post) / CZ < 0.0){
8708 std::cout <<
"DSQRT 2 < 0 " << CZ << std::endl;
8709 std::cout <<
"This event was not considered" << std::endl;
8713 ZA1width = std::sqrt(0.5*std::sqrt(1.0/A_levdens*E_scission_post)/CZ);
8724 ZA1width =
max(ZA1width,sigZmin);
8726 if(imode == 1 && cpol1 != 0.0){
8730 Z1rr = Z1UCD - cpol1 * A_scission/N_scission;
8736 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED" << std::endl;
8739 if ((
utilabs(Z1rr - Z1r) > 3.0*ZA1width) || Z1r<1.0)
goto fiss2801;
8742 if( imode == 2 && cpol2 != 0.0 ){
8746 Z1rr = Z1UCD - cpol2 * A_scission/N_scission;
8751 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED" << std::endl;
8754 if( (
utilabs(Z1rr - Z1r) > 3.0*ZA1width) || Z1r < 1.0 )
goto fiss2802;
8762 re1 =
frldm( Z1UCD-1.0, N1r+1.0, beta1 ) +
8763 frldm( Z2UCD+1.0, N2r-1.0, beta2 ) +
8764 ecoul( Z1UCD-1.0, N1r+1.0, beta1,
8765 Z2UCD+1.0, N2r-1.0, beta2, d );
8766 re2 =
frldm( Z1UCD, N1r, beta1) +
8767 frldm( Z2UCD, N2r, beta2 ) +
8768 ecoul( Z1UCD, N1r, beta1,
8769 Z2UCD, N2r, beta2, d );
8770 re3 =
frldm( Z1UCD+1.0, N1r-1.0, beta1 ) +
8771 frldm( Z2UCD-1.0, N2r+1.0, beta2 ) +
8772 ecoul( Z1UCD+1.0, N1r-1.0, beta1,
8773 Z2UCD-1.0, N2r+1.0, beta2, d );
8774 eps2 = ( re1 - 2.0*re2 + re3 ) / 2.0;
8775 eps1 = ( re3 - re1 ) / 2.0;
8776 DN1_POL = -eps1 / ( 2.0 * eps2 );
8778 Z1rr = Z1UCD + DN1_POL;
8783 DN1_POL = DN1_POL - 0.6 *
Uwash(E_scission_post,Ecrit,FREDSHELL,gamma);
8784 Z1rr = Z1UCD + DN1_POL;
8785 if ( Z1rr < 50. ) Z1rr = 50.0;
8787 DN1_POL = DN1_POL + 0.60 *
Uwash(E_scission_post,Ecrit,FREDSHELL,gamma);
8788 Z1rr = Z1UCD + DN1_POL;
8789 if ( Z1rr > 50.0 ) Z1rr = 50.0;
8799 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED" << std::endl;
8803 if( (
utilabs(Z1rr - Z1r) > 3.0*ZA1width) || (Z1r < 1.0) )
goto fiss2803;
8815 z2 =
dint( Z_scission ) - z1;
8817 N2 =
dint( N_scission ) - N1;
8821 if( (z1 < 0) || (z2 < 0) || (a1 < 0) || (a2 < 0) ){
8822 std::cout <<
" -------------------------------" << std::endl;
8823 std::cout <<
" Z, A, N : " <<
Z <<
" " <<
A <<
" " << N << std::endl;
8824 std::cout << z1 <<
" " << z2 <<
" " << a1 <<
" " << a2 << std::endl;
8825 std::cout << E_scission_post <<
" " << A_levdens <<
" " << CZ << std::endl;
8827 std::cout <<
" -------------------------------" << std::endl;
8836 if( N1mean > N*0.50 ){
8840 if(beta2< beta2gs) beta2 = beta2gs;
8841 E1exc = E_scission_pre * a1 /
A + E_defo;
8842 E_defo =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8843 E2exc = E_scission_pre * a2 /
A + E_defo;
8847 if(beta1< beta1gs) beta1 = beta1gs;
8848 E_defo =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8849 E1exc = E_scission_pre * a1 /
A + E_defo;
8851 E2exc = E_scission_pre * a2 /
A + E_defo;
8858 if( N1mean > N*0.5 ){
8861 if(beta1< beta1gs) beta1 = beta1gs;
8862 E_defo =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8863 E1exc = E_scission_pre * a1 /
A + E_defo;
8865 if(beta2< beta2gs) beta2 = beta2gs;
8866 E_defo =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8867 E2exc = E_scission_pre * a2 /
A + E_defo;
8871 if(beta2< beta2gs) beta2 = beta2gs;
8872 E_defo =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8873 E2exc = E_scission_pre * a2 /
A + E_defo;
8875 if(beta1< beta1gs) beta1 = beta1gs;
8876 E_defo =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8877 E1exc = E_scission_pre * a1 /
A + E_defo;
8884 if(beta1< beta1gs) beta1 = beta1gs;
8886 if(beta2< beta2gs) beta2 = beta2gs;
8887 E_defo1 =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8888 E_defo2 =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8889 E1exc = E_scission_pre * a1 /
A + E_defo1;
8890 E2exc = E_scission_pre * a2 /
A + E_defo2;
8895 TKER = ( z1 * z2 * 1.440 ) /
8896 ( R0 * std::pow(a1,0.333330) * (1.0 + 2.0/3.0 * beta1 ) +
8897 R0 * std::pow(a2,0.333330) * (1.0 + 2.0/3.0 * beta2 ) + 2.0 );
8899 EkinR1 = TKER * a2 /
A;
8900 EkinR2 = TKER * a1 /
A;
8901 v1 = std::sqrt(EkinR1/a1) * 1.3887;
8902 v2 = std::sqrt(EkinR2/a2) * 1.3887;
8911 if(
e1<0.)
goto fis987;
8915 if(
e2<0.)
goto fis988;
8917 (*NbLam0_par) = NbLam0;
8955 G4double r_in = 0.0, r_rest = 0.0, r_help = 0.0;
8961 r_in = r_origin + 0.5;
8963 if (r_even_odd < 0.001) {
8964 i_out = (
G4int)(r_floor);
8967 r_rest = r_in - r_floor;
8968 r_middle = r_floor + 0.5;
8969 n_floor = (
G4int)(r_floor);
8970 if (n_floor%2 == 0) {
8972 r_help = r_middle + (r_rest - 0.5) * (1.0 - r_even_odd);
8976 r_help = r_middle + (r_rest - 0.5) * (1.0 + r_even_odd);
8978 i_out = (
G4int)(r_help);
8994 G4double xcom = 0.0, xvs = 0.0, xe = 0.0;
9000 xcom = 1.0 - 1.7826 * ((a - 2.0*z)/a)*((a - 2.0*z)/a);
9002 xvs = - xcom * ( 15.4941 * a -
9003 17.9439 * std::pow(a,2.0/3.0) * (1.0+0.4*
alpha*
alpha) );
9005 xe = z*z * (0.7053/(std::pow(a,1.0/3.0)) * (1.0-0.2*
alpha*
alpha) - 1.1529/a);
9032 dtot = r0 * ( std::pow((z1+n1),1.0/3.0) * (1.0+0.6666667*beta1)
9033 + std::pow((z2+n2),1.0/3.0) * (1.0+0.6666667*beta2) ) + d;
9034 fecoul = z1 * z2 * 1.44 / dtot;
9046 R_wash = std::exp(-E * Freduction * gamma);
9048 R_wash = std::exp(- Ecrit * Freduction * gamma -(E-Ecrit) * gamma);
9106 G4double z = 0.0,
n = 0.0, a = 0.0, av = 0.0, as = 0.0;
9107 G4double a0 = 0.0, c1 = 0.0, c4 = 0.0, b1 = 0.0, b3 = 0.0;
9108 G4double ff = 0.0, ca = 0.0, w = 0.0, efl = 0.0;
9109 G4double r0 = 0.0, kf = 0.0, ks = 0.0;
9110 G4double kv = 0.0, rp = 0.0, ay = 0.0, aden = 0.0, x0 = 0.0, y0 = 0.0;
9111 G4double esq = 0.0, ael = 0.0, i = 0.0;
9162 c1 = 3.0/5.0*esq/r0;
9163 c4 = 5.0/4.0*std::pow((3.0/(2.0*
pi)),(2.0/3.0)) * c1;
9164 kf = std::pow((9.0*
pi*z/(4.0*a)),(1.0/3.0))/r0;
9166 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));
9170 x0 = r0 * std::pow(a,(1.0/3.0)) / ay;
9171 y0 = r0 * std::pow(a,(1.0/3.0)) / aden;
9173 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);
9175 b3 = 1.0 - 5.0/std::pow(y0,2) * (1.0 - 15.0/(8.0*y0) + 21.0/(8.0 * std::pow(y0,3))
9176 - 3.0/4.0 * (1.0 + 9.0/(2.0*y0) + 7.0/std::pow(y0,2)
9177 + 7.0/(2.0 * std::pow(y0,3))) * std::exp(-2.0*y0));
9181 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
9182 + 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))
9183 + ff*std::pow(z,2)/a -ca*(
n-z) - ael * std::pow(z,(2.39e0));
9189 return eflmacResult;
9194void 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){
9196 G4int INMIN,INMAX,NDIF=0,IMEM;
9197 G4int NEVA=0,PEVA=0;
9205 for(
G4int i=0;i<200;i++){
9206 BU_TAB_TEMP[i][0] = 0.0;
9207 BU_TAB_TEMP[i][1] = 0.0;
9208 BU_TAB_TEMP[i][2] = 0.0;
9209 BU_TAB_TEMP[i][3] = 0.0;
9210 BU_TAB_TEMP[i][4] = 0.0;
9217 if(AFP==0 && ZFP==0){
9221 if((AFP==1 && ZFP==0) ||
9222 (AFP==1 && ZFP==1) ||
9223 (AFP==2 && ZFP==1) ||
9224 (AFP==3 && ZFP==1) ||
9225 (AFP==3 && ZFP==2) ||
9226 (AFP==4 && ZFP==2) ||
9227 (AFP==6 && ZFP==2) ||
9236 if ((AFP-ZFP)==0 && ZFP>1){
9237 for(
G4int I = 0;I<=AFP-2;I++){
9239 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9240 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9241 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9242 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9243 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9244 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9245 *ILOOP = *ILOOP + 1;
9262 for(
G4int I = 1;I<=10; I++){
9274 for(
G4int I = 0;I< IMEM;I++){
9280 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9281 BU_TAB_TEMP[I+1+*ILOOP][0] = 1.0;
9282 BU_TAB_TEMP[I+1+*ILOOP][1] = 1.0;
9283 BU_TAB_TEMP[I+1+*ILOOP][2] = VP2X;
9284 BU_TAB_TEMP[I+1+*ILOOP][3] = VP2Y;
9285 BU_TAB_TEMP[I+1+*ILOOP][4] = VP2Z;
9290 *ILOOP = *ILOOP + IMEM;
9295 NEVA = NDIF - INMAX;
9298 for(
G4int I = 0;I<NEVA;I++){
9304 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9306 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9307 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9308 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9309 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9310 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9311 *ILOOP = *ILOOP + 1;
9318 if ((AFP>=2) && (ZFP==0)){
9319 for(
G4int I = 0;I<= AFP-2;I++){
9323 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9325 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9326 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9327 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9328 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9329 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9330 *ILOOP = *ILOOP + 1;
9342 std::cout <<
"WARNING - BU UNSTABLE: AF < ZF" << std::endl;
9347 if ((AFP>=4) && (ZFP==1)){
9349 for(
G4int I = 0; I<AFP-3;I++){
9353 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9355 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9356 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9357 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9358 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9359 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9360 *ILOOP = *ILOOP + 1;
9372 if ((AFP==4) && (ZFP==3)){
9380 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9382 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9383 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9384 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9385 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9386 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9387 *ILOOP = *ILOOP + 1;
9389 if ((AFP==5) && (ZFP==2)){
9397 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9398 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9399 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9400 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9401 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9402 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9403 *ILOOP = *ILOOP + 1;
9406 if ((AFP==5) && (ZFP==3)){
9414 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9415 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9416 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9417 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9418 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9419 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9420 *ILOOP = *ILOOP + 1;
9423 if ((AFP==6) && (ZFP==4)){
9432 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9433 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9434 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9435 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9436 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9437 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9438 *ILOOP = *ILOOP + 1;
9446 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9447 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9448 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9449 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9450 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9451 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9452 *ILOOP = *ILOOP + 1;
9454 if ((AFP==7)&&(ZFP==2)){
9462 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9463 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9464 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9465 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9466 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9467 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9468 *ILOOP = *ILOOP + 1;
9471 if ((AFP==7) && (ZFP==5)){
9473 for(
int I = 0; I<= AFP-5;I++){
9475 double(AFP-I-1),
double(ZFP-I-1),
9477 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9478 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9479 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9480 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9481 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9482 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9483 *ILOOP = *ILOOP + 1;
9494 if ((AFP==8) && (ZFP==4)){
9501 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9502 BU_TAB_TEMP[*ILOOP][0] = 2.0;
9503 BU_TAB_TEMP[*ILOOP][1] = 4.0;
9504 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9505 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9506 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9507 *ILOOP = *ILOOP + 1;
9509 if ((AFP==8) && (ZFP==6)){
9518 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9519 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9520 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9521 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9522 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9523 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9524 *ILOOP = *ILOOP + 1;
9531 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9532 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9533 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9534 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9535 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9536 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9537 *ILOOP = *ILOOP + 1;
9543 if((AFP==9) && (ZFP==2)){
9552 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9553 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9554 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9555 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9556 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9557 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9558 *ILOOP = *ILOOP + 1;
9564 if((AFP==9) && (ZFP==5)){
9572 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9573 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9574 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9575 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9576 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9577 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9578 *ILOOP = *ILOOP + 1;
9585 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9586 BU_TAB_TEMP[*ILOOP][0] = 2.0;
9587 BU_TAB_TEMP[*ILOOP][1] = 4.0;
9588 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9589 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9590 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9591 *ILOOP = *ILOOP + 1;
9597 if((AFP==10) && (ZFP==2)){
9606 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9607 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9608 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9609 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9610 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9611 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9612 *ILOOP = *ILOOP + 1;
9620 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9621 BU_TAB_TEMP[*ILOOP][0] = 0.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;
9631 if ((AFP==10) && (ZFP==3)){
9639 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9640 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9641 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9642 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9643 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9644 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9645 *ILOOP = *ILOOP + 1;
9650 if ((AFP==10) && (ZFP==7)){
9658 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9659 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9660 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9661 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9662 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9663 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9664 *ILOOP = *ILOOP + 1;
9670 if((AFP==11) && (ZFP==7)){
9678 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9679 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9680 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9681 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9682 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9683 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9684 *ILOOP = *ILOOP + 1;
9689 if ((AFP==12) && (ZFP==8)){
9698 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9699 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9700 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9701 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9702 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9703 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9704 *ILOOP = *ILOOP + 1;
9711 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9712 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9713 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9714 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9715 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9716 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9717 *ILOOP = *ILOOP + 1;
9722 if ((AFP==15) && (ZFP==9)){
9730 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9731 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9732 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9733 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9734 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9735 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9736 *ILOOP = *ILOOP + 1;
9742 if ((AFP==16) && (ZFP==9)){
9750 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9751 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9752 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9753 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9754 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9755 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9756 *ILOOP = *ILOOP + 1;
9762 if ((AFP==16) && (ZFP==10)){
9770 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9771 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9772 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9773 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9774 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9775 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9776 *ILOOP = *ILOOP + 1;
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((AFP==18) && (ZFP==11)){
9802 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9803 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9804 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9805 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9806 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9807 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9808 *ILOOP = *ILOOP + 1;
9813 if((AFP==19) && (ZFP==11)){
9821 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9822 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9823 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9824 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9825 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9826 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9827 *ILOOP = *ILOOP + 1;
9832 if (ZFP>=4 && (AFP-ZFP)==1){
9837 for(
G4int I = 0;I< NEVA;I++){
9841 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9842 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9843 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9844 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9845 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9846 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9847 *ILOOP = *ILOOP + 1;
9852 for(
G4int I = 0;I<PEVA;I++){
9856 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9857 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9858 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9859 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9860 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9861 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9862 *ILOOP = *ILOOP + 1;
9880void 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){
9883 G4double PX1,PX2,PY1,PY2,PZ1,PZ2,PTOT;
9884 G4double RNDT,CTET1,STET1,RNDP,PHI1,ETOT_P1,ETOT_P2;
9886 G4double vxout=0.,vyout=0.,vzout=0.;
9887 G4int iain,izin,ianew,iznew,inin,innew;
9897 innew = ianew - iznew;
9902 mglms(ain,zin,3,&MASS);
9903 mglms(anew,znew,3,&MASS1);
9904 mglms(ain-anew,zin-znew,3,&MASS2);
9905 ekin_tot = MASS-MASS1-MASS2;
9909 if(izin>12)std::cout <<
"*** ZIN > 12 ***" << izin << std::endl;
9912 if( ekin_tot<0.00 ){
9923 EKIN_P1 = ekin_tot*(ain-anew)/ ain;
9924 ETOT_P1 = EKIN_P1 + anew * AMU;
9925 PTOT = anew*AMU*std::sqrt((EKIN_P1/(anew*AMU)+1.0)*(EKIN_P1/(anew*AMU)+1.0)-1.0);
9928 CTET1 = 2.0*RNDT-1.0;
9929 STET1 = std::sqrt(1.0-CTET1*CTET1);
9931 PHI1 = RNDP*2.0*3.141592654;
9932 PX1 = PTOT * STET1*std::cos(PHI1);
9933 PY1 = PTOT * STET1*std::sin(PHI1);
9935 *v1x =
C * PX1 / ETOT_P1;
9936 *v1y =
C * PY1 / ETOT_P1;
9937 *v1z =
C * PZ1 / ETOT_P1;
9938 lorentz_boost(vxin,vyin,vzin,*v1x,*v1y,*v1z,&vxout,&vyout,&vzout);
9946 ETOT_P2 = (ekin_tot - EKIN_P1) + (ain-anew) * AMU;
9947 *v2x =
C * PX2 / ETOT_P2;
9948 *v2y =
C * PY2 / ETOT_P2;
9949 *v2z =
C * PZ2 / ETOT_P2;
9950 lorentz_boost(vxin,vyin,vzin,*v2x,*v2y,*v2z,&vxout,&vyout,&vzout);
9968 G4double GAMMA,VR,
C,CC,DENO,VXNOM,VYNOM,VZNOM;
9979 VR = std::sqrt(VXR*VXR + VYR*VYR + VZR*VZR);
9986 GAMMA = 1.0/std::sqrt(1.0 - VR*VR/CC);
9987 DENO = 1.0 - VXR*VXIN/CC - VYR*VYIN/CC - VZR*VZIN/CC;
9990 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;
9992 *VXOUT = VXNOM / (GAMMA * DENO);
9995 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;
9997 *VYOUT = VYNOM / (GAMMA * DENO);
10000 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;
10002 *VZOUT = VZNOM / (GAMMA * DENO);
10014 G4double EFF1=0.,EFF2=0.,VFF1=0.,VFF2=0.,
10015 AF1=0.,ZF1=0.,AF2=0.,ZF2=0.,
10016 AFF1=0.,ZFF1=0.,AFF2=0.,ZFF2=0.,
10017 vz1_eva=0., vx1_eva=0.,vy1_eva=0.,
10018 vz2_eva=0., vx2_eva=0.,vy2_eva=0.,
10019 vx_eva_sc=0.,vy_eva_sc=0.,vz_eva_sc=0.,
10020 VXOUT=0.,VYOUT=0.,VZOUT=0.,
10021 VX2OUT=0.,VY2OUT=0.,VZ2OUT=0.;
10022 G4int IEV_TAB_FIS=0,IEV_TAB_TEMP=0;
10023 G4double EV_TEMP1[200][6], EV_TEMP2[200][6],mtota=0.;
10024 G4int inttype = 0,inum=0;
10027 G4int NbLam0= (*NbLam0_par);
10029 for(
G4int I1=0;I1<200;I1++)
10030 for(
G4int I2=0;I2<6;I2++){
10031 EV_TEMP[I1][I2] = 0.0;
10032 EV_TEMP1[I1][I2] = 0.0;
10033 EV_TEMP2[I1][I2] = 0.0;
10036 G4double et = EE - JPRF * JPRF * 197. * 197./(2.*0.4*931.*std::pow(AF,5.0/3.0)*1.16*1.16);
10038 fissionDistri(AF,ZF,et,AF1,ZF1,EFF1,VFF1,AF2,ZF2,EFF2,VFF2,
10039 vx_eva_sc,vy_eva_sc,vz_eva_sc,&NbLam0);
10044 G4double pbH = (AF1 - ZF1) / (AF1 - ZF1 + AF2 - ZF2);
10045 for(
G4int i=0;i<NbLam0;i++){
10065 G4double VPERP1 = std::sqrt(VFF1*VFF1 - VZ1_FISSION*VZ1_FISSION);
10067 G4double VX1_FISSION = VPERP1 * std::sin(ALPHA1);
10068 G4double VY1_FISSION = VPERP1 * std::cos(ALPHA1);
10069 G4double VX2_FISSION = - VX1_FISSION / VFF1 * VFF2;
10070 G4double VY2_FISSION = - VY1_FISSION / VFF1 * VFF2;
10071 G4double VZ2_FISSION = - VZ1_FISSION / VFF1 * VFF2;
10074 if( (ZF1<=0.0) || (AF1<=0.0) || (AF1<ZF1) ){
10075 std::cout <<
"F1 unphysical: "<<ZF<<
" "<<AF<<
" "<<EE<<
" "<<ZF1<<
" "<<AF1 << std::endl;
10081 G4int FF11=0, FIMF11=0;
10082 G4double ZIMFF1=0., AIMFF1=0.,TKEIMF1=0.,JPRFOUT=0.;
10084 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);
10086 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
10087 EV_TEMP[IJ+IEV_TAB_FIS][0] = EV_TEMP1[IJ][0];
10088 EV_TEMP[IJ+IEV_TAB_FIS][1] = EV_TEMP1[IJ][1];
10095 EV_TEMP1[IJ][2],EV_TEMP1[IJ][3],EV_TEMP1[IJ][4],
10096 &VXOUT,&VYOUT,&VZOUT);
10099 &VX2OUT,&VY2OUT,&VZ2OUT);
10100 EV_TEMP[IJ+IEV_TAB_FIS][2] = VX2OUT;
10101 EV_TEMP[IJ+IEV_TAB_FIS][3] = VY2OUT;
10102 EV_TEMP[IJ+IEV_TAB_FIS][4] = VZ2OUT;
10105 IEV_TAB_FIS = IEV_TAB_FIS + IEV_TAB_TEMP;
10110 if( (ZF2<=0.0) || (AF2<=0.0) || (AF2<ZF2) ){
10111 std::cout <<
"F2 unphysical: "<<ZF<<
" "<<AF<<
" "<<EE<<
" "<<ZF2<<
" "<<AF2 << std::endl;
10117 G4int FF22=0, FIMF22=0;
10118 G4double ZIMFF2=0., AIMFF2=0.,TKEIMF2=0.,JPRFOUT=0.;
10120 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);
10122 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
10123 EV_TEMP[IJ+IEV_TAB_FIS][0] = EV_TEMP2[IJ][0];
10124 EV_TEMP[IJ+IEV_TAB_FIS][1] = EV_TEMP2[IJ][1];
10131 EV_TEMP2[IJ][2],EV_TEMP2[IJ][3],EV_TEMP2[IJ][4],
10132 &VXOUT,&VYOUT,&VZOUT);
10135 &VX2OUT,&VY2OUT,&VZ2OUT);
10136 EV_TEMP[IJ+IEV_TAB_FIS][2] = VX2OUT;
10137 EV_TEMP[IJ+IEV_TAB_FIS][3] = VY2OUT;
10138 EV_TEMP[IJ+IEV_TAB_FIS][4] = VZ2OUT;
10141 IEV_TAB_FIS = IEV_TAB_FIS + IEV_TAB_TEMP;
10154 VX1_FISSION,VY1_FISSION,VZ1_FISSION,
10155 &VXOUT,&VYOUT,&VZOUT);
10156 VX1_FISSION = VXOUT;
10157 VY1_FISSION = VYOUT;
10158 VZ1_FISSION = VZOUT;
10160 VX2_FISSION,VY2_FISSION,VZ2_FISSION,
10161 &VXOUT,&VYOUT,&VZOUT);
10162 VX2_FISSION = VXOUT;
10163 VY2_FISSION = VYOUT;
10164 VZ2_FISSION = VZOUT;
10169 (*VX1_FISSION_par) = VX1_FISSION;
10170 (*VY1_FISSION_par) = VY1_FISSION;
10171 (*VZ1_FISSION_par) = VZ1_FISSION;
10172 (*VX_EVA_SC_par)=vx_eva_sc;
10173 (*VY_EVA_SC_par)=vy_eva_sc;
10174 (*VZ_EVA_SC_par)=vz_eva_sc;
10178 (*VX2_FISSION_par) = VX2_FISSION;
10179 (*VY2_FISSION_par) = VY2_FISSION;
10180 (*VZ2_FISSION_par) = VZ2_FISSION;
10181 (*IEV_TAB_FIS_par) = IEV_TAB_FIS;
10182 (*NbLam0_par) = NbLam1 + NbLam2;
10183 if(NbLam0>(NbLam1 + NbLam2))
varntp->
kfis = 25;
10190 G4double V_over_V0,R0,RALL,RHAZ,R,TKE,Ekin,V,VPERP,ALPHA1;
10202 RALL = R0 * std::pow(V_over_V0,1.0/3.0) * std::pow(AAL,1.0/3.0);
10204 R = std::pow(RHAZ,1.0/3.0) * RALL;
10205 TKE = 1.44 *
Z * ZALL * R*R * (1.0 -
A/AAL)*(1.0 -
A/AAL) / std::pow(RALL,3.0);
10207 Ekin = TKE * (AAL -
A) / AAL;
10209 V = std::sqrt(Ekin/
A) * 1.3887;
10211 VPERP = std::sqrt(V*V - (*VZ)*(*VZ));
10213 *VX = VPERP * std::sin(ALPHA1);
10214 *VY = VPERP * std::cos(ALPHA1);
10240 ix =
G4int(y * 100 + 43543000);
10241 if(
mod(ix,2) == 0) {
10259 y=std::pow(
G4AblaRandom::flat()*(std::pow(rxmax,l_plus)-std::pow(rxmin,l_plus))+ std::pow(rxmin,l_plus),1.0/l_plus);
10293 G4double Efermi = 5.0 * SIGMA_0 * SIGMA_0 / (2.0 * 931.4940);
10299 SIGMA_0 = SIGMA_0 * std::pow(V0_over_VBU,1.0/3.0);
10304 GOLDHA_BU = SIGMA_0 * std::sqrt((APRF*(AABRA-APRF))/(AABRA-1.0));
10305 GOLDHA = GOLDHA_BU*std::sqrt(1.0 +
10313 (AABRA - APRF) / AABRA);
10314 GOLDHA = GOLDHA_BU;
10319 GOLDHA = SIGMA_0 * std::sqrt((APRF*(AABRA-APRF))/(AABRA-1.0));
10327 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING PX IN Rn07.FOR. A VALUE WILL BE FORCED." << std::endl;
10328 *PX = (AABRA-1.0)*931.4940;
10330 if(std::abs(*PX)>= AABRA*931.494){
10339 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING PY IN Rn07.FOR. A VALUE WILL BE FORCED." << std::endl;
10340 *PY = (AABRA-1.0)*931.4940;
10342 if(std::abs(*PY)>= AABRA*931.494){
10351 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING PZ IN Rn07.FOR. A VALUE WILL BE FORCED." << std::endl;
10352 *PZ = (AABRA-1.0)*931.4940;
10354 if(std::abs(*PZ)>= AABRA*931.494){
10371 v1 = 2.0*
haz(k) - 1.0;
10372 v2 = 2.0*
haz(k) - 1.0;
10373 r = std::pow(v1,2) + std::pow(v2,2);
10376 fac = std::sqrt(-2.*std::log(r)/r);
10378 fgausshaz = v2*
fac*sig+xmoy;
10382 fgausshaz=gset*sig+xmoy;
static const G4double e1[44]
static const G4double e2[44]
G4double C(G4double temp)
G4double B(G4double temperature)
G4double D(G4double temp)
G4double Y(G4double density)
static const G4double eps
static const G4int amax[]
static const G4int amin[]
static const G4double alpha
static const G4double fac
static constexpr double m
static constexpr double pc
static constexpr double pi
G4double getMexp(G4int A, G4int Z)
G4double getPace2(G4int A, G4int Z)
G4double getAlpha(G4int A, G4int Z)
G4double getBeta2(G4int A, G4int Z)
G4double getRms(G4int A, G4int Z)
G4double getVgsld(G4int A, G4int Z)
G4double getBeta4(G4int A, G4int Z)
G4double getEcnz(G4int A, G4int Z)
G4int getMexpID(G4int A, G4int Z)
G4double eflmac(G4int ia, G4int iz, G4int flag, G4int optshp)
G4double func_trans(G4double TIME, G4double ZF, G4double AF, G4double BET, G4double Y, G4double FT, G4double T_0)
void part_fiss(G4double BET, G4double GP, G4double GF, G4double Y, G4double TAUF, G4double TS1, G4double TSUM, G4int *CHOICE, G4double ZF, G4double AF, G4double FT, G4double *T_LAPSE, G4double *GF_LOC)
G4int IPOWERLIMHAZ(G4double lambda, G4int xmin, G4int xmax)
void imf(G4double ACN, G4double ZCN, G4double TEMP, G4double EE, G4double *ZIMF, G4double *AIMF, G4double *BIMF, G4double *SBIMF, G4double *TIMF, G4double JPRF)
G4double fvmaxhaz_neut(G4double x)
G4double tau(G4double bet, G4double homega, G4double ef, G4double t)
G4double fmaxhaz(G4double T)
G4double gammp(G4double a, G4double x)
void isostab_lim(G4int z, G4int *nmin, G4int *nmax)
void barrs(G4int Z1, G4int A1, G4int Z2, G4int A2, G4double *sBARR, G4double *sOMEGA)
void even_odd(G4double r_origin, G4double r_even_odd, G4int &i_out)
G4double DSIGN(G4double a, G4double b)
G4double fvmaxhaz(G4double T)
G4double gethyperseparation(G4double A, G4double Z, G4int ny)
G4double EV_TAB_SSC[200][6]
G4int nint(G4double number)
G4double gammln(G4double xx)
void fomega_gs(G4double AF, G4double ZF, G4double *K1, G4double *sOMEGA, G4double *sHOMEGA)
G4double tunnelling(G4double A, G4double ZPRF, G4double Y, G4double EE, G4double EF, G4double TEMP, G4double DENSG, G4double DENSF, G4double ENH_FACT)
G4int ISIGN(G4int a, G4int b)
G4double fmaxhaz_old(G4double T)
void evap_postsaddle(G4double A, G4double Z, G4double E_scission_pre, G4double *E_scission_post, G4double *A_scission, G4double *Z_scission, G4double &vx_eva, G4double &vy_eva, G4double &vz_eva, G4int *NbLam0_par)
G4double pace2(G4double a, G4double z)
void mglw(G4double a, G4double z, G4double *el)
G4double dmin1(G4double a, G4double b, G4double c)
G4int mod(G4int a, G4int b)
void unstable_tke(G4double AIN, G4double ZIN, G4double ANEW, G4double ZNEW, G4double VXIN, G4double VYIN, G4double VZIN, G4double *V1X, G4double *V1Y, G4double *V1Z, G4double *V2X, G4double *V2Y, G4double *V2Z)
void parite(G4double n, G4double *par)
G4Abla(G4Volant *aVolant, G4VarNtp *aVarntp)
G4double fissility(G4int a, G4int z, G4int ny, G4double sn, G4double slam, G4int optxfis)
void unbound(G4double SN, G4double SP, G4double SD, G4double ST, G4double SHE, G4double SA, G4double BP, G4double BD, G4double BT, G4double BHE, G4double BA, G4double *PROBF, G4double *PROBN, G4double *PROBP, G4double *PROBD, G4double *PROBT, G4double *PROBHE, G4double *PROBA, G4double *PROBIMF, G4double *PROBG, G4double *ECN, G4double *ECP, G4double *ECD, G4double *ECT, G4double *ECHE, G4double *ECA)
G4double cram(G4double bet, G4double homega)
void AMOMENT(G4double AABRA, G4double APRF, G4int IMULTIFR, G4double *PX, G4double *PY, G4double *PZ)
void unstable_nuclei(G4int AFP, G4int ZFP, G4int *AFPNEW, G4int *ZFPNEW, G4int &IOUNSTABLE, G4double VX, G4double VY, G4double VZ, G4double *VP1X, G4double *VP1Y, G4double *VP1Z, G4double BU_TAB_TEMP[200][6], G4int *ILOOP)
void tke_bu(G4double Z, G4double A, G4double ZALL, G4double AAL, G4double *VX, G4double *VY, G4double *VZ)
void mglms(G4double a, G4double z, G4int refopt4, G4double *el)
void barfit(G4int iz, G4int ia, G4int il, G4double *sbfis, G4double *segs, G4double *selmax)
G4double dint(G4double a)
G4double umass(G4double z, G4double n, G4double beta)
G4int idnint(G4double value)
void appariem(G4double a, G4double z, G4double *del)
void gser(G4double *gamser, G4double a, G4double x, G4double gln)
G4double spdef(G4int a, G4int z, G4int optxfis)
void qrot(G4double z, G4double a, G4double bet, G4double sig, G4double u, G4double *qr)
G4double utilabs(G4double a)
void direct(G4double zprf, G4double a, G4double ee, G4double jprf, G4double *probp_par, G4double *probd_par, G4double *probt_par, G4double *probn_par, G4double *probhe_par, G4double *proba_par, G4double *probg_par, G4double *probimf_par, G4double *probf_par, G4double *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)
void lpoly(G4double x, G4int n, G4double pl[])
void fission_width(G4double ZPRF, G4double A, G4double EE, G4double BS, G4double BK, G4double EF, G4double Y, G4double *GF, G4double *TEMP, G4double JPR, G4int IEROT, G4int FF_ALLOWED, G4int OPTCOL, G4int OPTSHP, G4double DENSG)
G4double Uwash(G4double E, G4double Ecrit, G4double Freduction, G4double gamma)
G4int max(G4int a, G4int b)
G4int min(G4int a, G4int b)
void fissionDistri(G4double &a, G4double &z, G4double &e, G4double &a1, G4double &z1, G4double &e1, G4double &v1, G4double &a2, G4double &z2, G4double &e2, G4double &v2, G4double &vx_eva_sc, G4double &vy_eva_sc, G4double &vz_eva_sc, G4int *NbLam0_par)
G4double width(G4double AMOTHER, G4double ZMOTHER, G4double APART, G4double ZPART, G4double TEMP, G4double B1, G4double SB1, G4double EXC)
G4double gethyperbinding(G4double A, G4double Z, G4int ny)
void guet(G4double *x_par, G4double *z_par, G4double *find_par)
G4double eflmac_profi(G4double a, G4double z)
G4double pen(G4double A, G4double ap, G4double omega, G4double T)
G4double ecoul(G4double z1, G4double n1, G4double beta1, G4double z2, G4double n2, G4double beta2, G4double d)
G4double gausshaz(G4int k, G4double xmoy, G4double sig)
G4double bipol(G4int iflag, G4double y)
void DeexcitationAblaxx(G4int nucleusA, G4int nucleusZ, G4double excitationEnergy, G4double angularMomentum, G4double momX, G4double momY, G4double momZ, G4int eventnumber)
void densniv(G4double a, G4double z, G4double ee, G4double ef, G4double *dens, G4double bshell, G4double bs, G4double bk, G4double *temp, G4int optshp, G4int optcol, G4double defbet, G4double *ecor, G4double jprf, G4int ifis, G4double *qr)
void fomega_sp(G4double AF, G4double Y, G4double *MFCD, G4double *sOMEGA, G4double *sHOMEGA)
void evapora(G4double zprf, G4double aprf, G4double *ee_par, G4double jprf, G4double *zf_par, G4double *af_par, G4double *mtota_par, G4double *vleva_par, G4double *vxeva_par, G4double *vyeva_par, G4int *ff_par, G4int *fimf_par, G4double *fzimf, G4double *faimf, G4double *tkeimf_par, G4double *jprfout, G4int *inttype_par, G4int *inum_par, G4double EV_TEMP[200][6], G4int *iev_tab_temp_par, G4int *nblam0)
G4double getdeltabinding(G4double a, G4int nblamb)
void FillData(G4int IMULTBU, G4int IEV_TAB)
void SetParametersG4(G4int z, G4int a)
void bsbkbc(G4double A, G4double Z, G4double *BS, G4double *BK, G4double *BC)
void fission(G4double AF, G4double ZF, G4double EE, G4double JPRF, G4double *VX1_FISSION, G4double *VY1_FISSION, G4double *VZ1_FISSION, G4double *VX2_FISSION, G4double *VY2_FISSION, G4double *VZ2_FISSION, G4int *ZFP1, G4int *AFP1, G4int *SFP1, G4int *ZFP2, G4int *AFP2, G4int *SFP2, G4int *imode, G4double *VX_EVA_SC, G4double *VY_EVA_SC, G4double *VZ_EVA_SC, G4double EV_TEMP[200][6], G4int *IEV_TAB_FIS, G4int *NbLam0)
void gcf(G4double *gammcf, G4double a, G4double x, G4double gln)
void lorentz_boost(G4double VXRIN, G4double VYRIN, G4double VZRIN, G4double VXIN, G4double VYIN, G4double VZIN, G4double *VXOUT, G4double *VYOUT, G4double *VZOUT)
void lorb(G4double AMOTHER, G4double ADAUGHTER, G4double LMOTHER, G4double EEFINAL, G4double *LORBITAL, G4double *SIGMA_LORBITAL)
void setVerboseLevel(G4int level)
G4double frldm(G4double z, G4double n, G4double beta)
G4double ecnz[EC2SUBROWS][EC2SUBCOLS]
G4double vgsld[ECLDROWS][ECLDCOLS]
G4double ecfnz[ECLDROWS][ECLDCOLS]
G4double alpha[ECLDROWS][ECLDCOLS]
G4double ecgnz[ECLDROWS][ECLDCOLS]
G4double rms[ECLDROWS][ECLDCOLS]
G4double beta2[ECLDROWSbeta][ECLDCOLSbeta]
G4double beta4[ECLDROWSbeta][ECLDCOLSbeta]
G4double efa[FBCOLS][FBROWS]
G4double bind[MASSIZEROWS][MASSIZECOLS]
G4double massexp[MASSIZEROWS][MASSIZECOLS]
G4int mexpiop[MASSIZEROWS][MASSIZECOLS]
G4double dm[PACESIZEROWS][PACESIZECOLS]
G4double enerj[VARNTPSIZE]
G4double pylab[VARNTPSIZE]
G4double pzlab[VARNTPSIZE]
G4int itypcasc[VARNTPSIZE]
G4double pxlab[VARNTPSIZE]
static const G4double Z1[5]