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++)
188 BU_TAB[I1][I2] = 0.0;
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;
193 EV_TAB[I1][I2] = 0.0;
194 EV_TAB_SSC[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;
256 if(T_freeze_out_in >= 0.0){
257 T_freeze_out = T_freeze_out_in;
259 T_freeze_out =
max(9.33*std::exp(-0.00282*AAINCL),5.5);
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);
269 T_diff = T_init - T_freeze_out;
271 if(T_diff>0.1 && zprf>2. && (aprf-zprf)>0.){
276 for(
G4int i=0;i<5;i++){
277 EE_diff = EINCL - a_tilda * T_freeze_out*T_freeze_out;
283 A_diff =
dint(EE_diff / (8.0 * 5.0 / T_freeze_out));
285 if(A_diff>AAINCL) A_diff = AAINCL;
287 A_FINAL = AAINCL - A_diff;
289 a_tilda = ald->
av*A_FINAL + ald->
as*std::pow(A_FINAL,2.0/3.0) + ald->
ak*std::pow(A_FINAL,1.0/3.0);
290 E_FINAL = a_tilda * T_freeze_out*T_freeze_out;
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);
332 E_FINAL = a_tilda * T_freeze_out*T_freeze_out;
334 ABU_SLOPE = (ASLOPE1-ASLOPE2)/4.0*(E_FINAL/AAINCL)+
335 ASLOPE1-(ASLOPE1-ASLOPE2)*7.0/4.0;
347 if(ABU_SLOPE > -1.01) ABU_SLOPE = -1.01;
350 Z_Breakup_sum = Z_FINAL;
354 for(
G4int i=0;i<100;i++){
361 std::cout <<
"WARNING: IPOWERLIMHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING A_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED: " << A_Breakup << std::endl;
365 if(A_Breakup>AAINCL)
goto mult4326;
368 std::cout <<
"A_BREAKUP <= 0 " << std::endl;
372 A_ACC = A_ACC + A_Breakup;
376 Z_Breakup_Mean =
dint(A_Breakup * ZAINCL / AAINCL);
378 Z_Breakup_sum = Z_Breakup_sum + Z_Breakup_Mean;
381 G_SYMM = 34.2281 - 5.14037 * E_FINAL/AAINCL;
382 if(E_FINAL/AAINCL < 2.0) G_SYMM = 25.0;
383 if(E_FINAL/AAINCL > 4.0) G_SYMM = 15.0;
388 CZ = 2.0 * G_SYMM * 4.0 / A_Breakup;
392 Sigma_Z = std::sqrt(T_freeze_out/CZ);
400 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING Z_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED: " << A_Breakup <<
" " << Z_Breakup << std::endl;
404 if(Z_Breakup<0.0 )
goto mult4333;
405 if((A_Breakup-Z_Breakup)<0.0)
goto mult4333;
406 if((A_Breakup-Z_Breakup)==0.0 && Z_Breakup!=1.0)
goto mult4333;
408 if(Z_Breakup>=ZAINCL){
411 std::cout <<
"Z_BREAKUP RESAMPLED MORE THAN 10 TIMES; EVENT WILL BE RESAMPLED AGAIN " << std::endl;
421 if(
idnint(A_Breakup-Z_Breakup)<INMIN ||
idnint(A_Breakup-Z_Breakup)>(INMAX+5)){
433 N_Breakup = A_Breakup - Z_Breakup;
434 BU_TAB[I_Breakup][0] =
dint(Z_Breakup);
435 BU_TAB[I_Breakup][1] =
dint(A_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){
520 BU_TAB[IEL][0] =
dint(ZTEMP);
521 BU_TAB[IEL][1] =
dint(ZTEMP + N_Breakup);
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;
551 BU_TAB[IEL][1] =
dint(BU_TAB[IEL][1]-
G4double(NBU_DIFF));
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){
569 ATEMP =
dint(BU_TAB[IEL][0] + N_Breakup);
570 }
else if(IEL>IMULTBU){
571 ATEMP =
dint(zprf + N_Breakup);
573 if((ATEMP - N_Breakup)<1.0 && N_Breakup<1.0){
583 BU_TAB[IEL][1] =
dint(BU_TAB[IEL][0] + N_Breakup);
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;
612 N_Breakup =
dint(BU_TAB[IEL][1]-BU_TAB[IEL][0]);
613 BU_TAB[IEL][0] =
dint(BU_TAB[IEL][0] -
G4double(ZBU_DIFF));
614 BU_TAB[IEL][1] =
dint(BU_TAB[IEL][0] + N_Breakup);
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;
627 N_Breakup =
dint(BU_TAB[IEL][1]-BU_TAB[IEL][0]);
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){
647 BU_TAB[IEL][0] =
dint(ZTEMP);
648 BU_TAB[IEL][1] =
dint(ZTEMP + N_Breakup);
652 aprf =
dint(ZTEMP + N_Breakup);
655 ZBU_DIFF = ZBU_DIFF -
ISIGN(1,ZBU_DIFF);
665 for(
G4int i =0;i<IMULTBU;i++){
670 if(BU_TAB[i][0]>2.0){
671 a_tilda_BU = ald->
av*BU_TAB[i][1] + ald->
as*std::pow(BU_TAB[i][1],2.0/3.0) + ald->
ak*std::pow(BU_TAB[i][1],1.0/3.0);
672 BU_TAB[i][2] = a_tilda_BU * T_freeze_out*T_freeze_out;
677 if(BU_TAB[i][0] > ZMEM){
687 BU_TAB[IMEM][0] = zprf;
688 BU_TAB[IMEM][1] = aprf;
689 BU_TAB[IMEM][2] = ee;
690 BU_TAB[IMEM][3] = jprf;
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++){
765 Bvalue_BU = Bvalue_BU +
eflmac(
idnint(BU_TAB[I_Breakup][1]),
idnint(BU_TAB[I_Breakup][0]),1,0);
766 Eexc_BU_SUM = Eexc_BU_SUM + BU_TAB[I_Breakup][2];
768 AMOMENT(AAINCL,BU_TAB[I_Breakup][1],1,&PX_BU,&PY_BU,&PZ_BU);
769 P_BU = std::sqrt(PX_BU*PX_BU + PY_BU*PY_BU + PZ_BU*PZ_BU);
772 ETOT_BU = std::sqrt(P_BU*P_BU + (BU_TAB[I_Breakup][1]*amu)*(BU_TAB[I_Breakup][1]*amu));
773 BU_TAB[I_Breakup][4] =
C * PX_BU / ETOT_BU;
774 BU_TAB[I_Breakup][5] =
C * PY_BU / ETOT_BU;
775 BU_TAB[I_Breakup][6] =
C * PZ_BU / ETOT_BU;
777 tke_bu(BU_TAB[I_Breakup][0],BU_TAB[I_Breakup][1],ZAINCL,AAINCL,&VX2_BU,&VY2_BU,&VZ2_BU);
784 BU_TAB[I_Breakup][4],BU_TAB[I_Breakup][5],BU_TAB[I_Breakup][6],
785 &VXOUT,&VYOUT,&VZOUT);
787 BU_TAB[I_Breakup][4] = VXOUT;
788 BU_TAB[I_Breakup][5] = VYOUT;
789 BU_TAB[I_Breakup][6] = VZOUT;
792 VBU2 = BU_TAB[I_Breakup][4]*BU_TAB[I_Breakup][4] +
793 BU_TAB[I_Breakup][5]*BU_TAB[I_Breakup][5] +
794 BU_TAB[I_Breakup][6]*BU_TAB[I_Breakup][6];
795 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C*
C));
796 ETOT_BU = BU_TAB[I_Breakup][1]*amu/GAMMA_REL;
797 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
798 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
799 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
801 PX_BU_SUM = PX_BU_SUM + PX_BU;
802 PY_BU_SUM = PY_BU_SUM + PY_BU;
803 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
808 P_BU_SUM = std::sqrt(PX_BU_SUM*PX_BU_SUM + PY_BU_SUM*PY_BU_SUM +
809 PZ_BU_SUM*PZ_BU_SUM);
812 ETOT_SUM = std::sqrt(P_BU_SUM*P_BU_SUM +
813 (AAINCL * amu)*(AAINCL * amu));
815 VX_BU_SUM =
C * PX_BU_SUM / ETOT_SUM;
816 VY_BU_SUM =
C * PY_BU_SUM / ETOT_SUM;
817 VZ_BU_SUM =
C * PZ_BU_SUM / ETOT_SUM;
825 VX_PREF,VY_PREF,VZ_PREF,
826 &VXOUT,&VYOUT,&VZOUT);
832 VBU2 = VX_PREF*VX_PREF + VY_PREF*VY_PREF + VZ_PREF*VZ_PREF;
833 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C*
C));
834 ETOT_PRF = aprf * amu / GAMMA_REL;
835 PXPRFP = ETOT_PRF * VX_PREF /
C;
836 PYPRFP = ETOT_PRF * VY_PREF /
C;
837 PZPRFP = ETOT_PRF * VZ_PREF /
C;
848 EKIN_BU = aprf * amu / GAMMA_REL - aprf * amu;
850 for(I_Breakup=0;I_Breakup<IMULTBU;I_Breakup++){
857 BU_TAB[I_Breakup][4],BU_TAB[I_Breakup][5],BU_TAB[I_Breakup][6],
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;
864 VBU2 = BU_TAB[I_Breakup][4]*BU_TAB[I_Breakup][4] +
865 BU_TAB[I_Breakup][5]*BU_TAB[I_Breakup][5] +
866 BU_TAB[I_Breakup][6]*BU_TAB[I_Breakup][6];
867 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C*
C));
869 ETOT_BU = BU_TAB[I_Breakup][1]*amu/GAMMA_REL;
871 EKIN_BU = EKIN_BU + BU_TAB[I_Breakup][1] * amu /
872 GAMMA_REL - BU_TAB[I_Breakup][1] * amu;
874 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
875 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
876 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
877 E_tot_BU = E_tot_BU + ETOT_BU;
879 PX_BU_SUM = PX_BU_SUM + PX_BU;
880 PY_BU_SUM = PY_BU_SUM + PY_BU;
881 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
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++){
937 BU_TAB[I_Breakup][4],BU_TAB[I_Breakup][5],BU_TAB[I_Breakup][6],
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;
944 VBU2 = BU_TAB[I_Breakup][4]*BU_TAB[I_Breakup][4] +
945 BU_TAB[I_Breakup][5]*BU_TAB[I_Breakup][5] +
946 BU_TAB[I_Breakup][6]*BU_TAB[I_Breakup][6];
947 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C*
C));
949 ETOT_BU = BU_TAB[I_Breakup][1]*amu/GAMMA_REL;
951 EKIN_BU = EKIN_BU + BU_TAB[I_Breakup][1] * amu /
952 GAMMA_REL - BU_TAB[I_Breakup][1] * amu;
954 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
955 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
956 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
957 E_tot_BU = E_tot_BU + ETOT_BU;
959 PX_BU_SUM = PX_BU_SUM + PX_BU;
960 PY_BU_SUM = PY_BU_SUM + PY_BU;
961 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
969 for(
G4int i=0;i<IMULTBU;i++){
970 if(BU_TAB[i][0]<3.0 || BU_TAB[i][0]==BU_TAB[i][1]){
972 BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6],
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;
1012 for(
G4int i=0;i<IMULTBU;i++)sumN=sumN+BU_TAB[i][1]-BU_TAB[i][0];
1014 for(
G4int i=0;i<IMULTBU;i++){
1015 problamb[i] = (BU_TAB[i][1]-BU_TAB[i][0])/sumN;
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++){
1041 EEBU = BU_TAB[i][2];
1042 BU_TAB[i][10] = BU_TAB[i][6];
1044 if(BU_TAB[i][0]>2.0){
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);
1049 BU_TAB[i][9] = jprfbu;
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;
1077 BU_TAB[i][4],BU_TAB[i][5],BU_TAB[i][6],
1078 &VXOUT,&VYOUT,&VZOUT);
1079 BU_TAB[i][4] = VXOUT;
1080 BU_TAB[i][5] = VYOUT;
1081 BU_TAB[i][6] = VZOUT;
1084 BU_TAB[i][7] =
dint(ZFBU);
1085 BU_TAB[i][8] =
dint(AFBU);
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;
1126 BU_TAB[i][4],BU_TAB[i][5],BU_TAB[i][6],
1127 &VXOUT,&VYOUT,&VZOUT);
1128 BU_TAB[i][4] = VXOUT;
1129 BU_TAB[i][5] = VYOUT;
1130 BU_TAB[i][6] = 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;
1166 BU_TAB[i][7] =
dint(ZFFBU);
1167 BU_TAB[i][8] =
dint(AFFBU);
1168 BU_TAB[i][11]= NbLamH;
1174 BU_TAB[i][4],BU_TAB[i][5],BU_TAB[i][6],
1175 &VXOUT,&VYOUT,&VZOUT);
1176 BU_TAB[i][4] = VXOUT;
1177 BU_TAB[i][5] = VYOUT;
1178 BU_TAB[i][6] = 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;
1210 BU_TAB[IMULTBU+ILOOPBU][0] = BU_TAB[i][0];
1211 BU_TAB[IMULTBU+ILOOPBU][1] = BU_TAB[i][1];
1212 BU_TAB[IMULTBU+ILOOPBU][2] = BU_TAB[i][2];
1213 BU_TAB[IMULTBU+ILOOPBU][3] = BU_TAB[i][3];
1214 BU_TAB[IMULTBU+ILOOPBU][7] =
dint(zffimf);
1215 BU_TAB[IMULTBU+ILOOPBU][8] =
dint(affimf);
1216 BU_TAB[IMULTBU+ILOOPBU][11]= NbLamimf;
1219 BU_TAB[i][4],BU_TAB[i][5],BU_TAB[i][6],
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;
1235 BU_TAB[i][7] = BU_TAB[i][0];
1236 BU_TAB[i][8] = BU_TAB[i][1];
1240 BU_TAB[i][11]= Nblamb[i];
1244 IMULTBU = IMULTBU + ILOOPBU;
1252 for(
G4int i=0;i<IMULTBU;i++){
1253 ABU_SUM = ABU_SUM + BU_TAB[i][8];
1254 ZBU_SUM = ZBU_SUM + BU_TAB[i][7];
1256 BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6],
1257 &VP1X,&VP1Y,&VP1Z,BU_TAB_TEMP1,&ILOOP);
1265 ABU_SUM = ABU_SUM +
G4double(afpnew) - BU_TAB[i][8];
1266 ZBU_SUM = ZBU_SUM +
G4double(zfpnew) - BU_TAB[i][7];
1269 BU_TAB[i][4] = VP1X;
1270 BU_TAB[i][5] = VP1Y;
1271 BU_TAB[i][6] = VP1Z;
1274 for(
G4int IJ=0;IJ<ILOOP;IJ++){
1275 BU_TAB[IMULTBU+INEWLOOP+IJ][7] = BU_TAB_TEMP1[IJ][0];
1276 BU_TAB[IMULTBU+INEWLOOP+IJ][8] = BU_TAB_TEMP1[IJ][1];
1277 BU_TAB[IMULTBU+INEWLOOP+IJ][4] = BU_TAB_TEMP1[IJ][2];
1278 BU_TAB[IMULTBU+INEWLOOP+IJ][5] = BU_TAB_TEMP1[IJ][3];
1279 BU_TAB[IMULTBU+INEWLOOP+IJ][6] = BU_TAB_TEMP1[IJ][4];
1280 BU_TAB[IMULTBU+INEWLOOP+IJ][2] = 0.0;
1281 BU_TAB[IMULTBU+INEWLOOP+IJ][3] = 0.0;
1282 BU_TAB[IMULTBU+INEWLOOP+IJ][0] = BU_TAB[i][0];
1283 BU_TAB[IMULTBU+INEWLOOP+IJ][1] = BU_TAB[i][1];
1284 BU_TAB[IMULTBU+INEWLOOP+IJ][11] = BU_TAB[i][11];
1285 ABU_SUM = ABU_SUM + BU_TAB[IMULTBU+INEWLOOP+IJ][8];
1286 ZBU_SUM = ZBU_SUM + BU_TAB[IMULTBU+INEWLOOP+IJ][7];
1289 INEWLOOP = INEWLOOP + ILOOP;
1294 IMULTBU = IMULTBU + INEWLOOP;
1298 VX_PREF,VY_PREF,VZ_PREF,
1299 &VXOUT,&VYOUT,&VZOUT);
1304 for(
G4int i=0;i<IMULTBU;i++){
1306 BU_TAB[i][4],BU_TAB[i][5],BU_TAB[i][6],
1307 &VXOUT,&VYOUT,&VZOUT);
1308 BU_TAB[i][4] = VXOUT;
1309 BU_TAB[i][5] = VYOUT;
1310 BU_TAB[i][6] = VZOUT;
1312 for(
G4int i=0;i<IEV_TAB;i++){
1314 EV_TAB[i][2],EV_TAB[i][3],EV_TAB[i][4],
1315 &VXOUT,&VYOUT,&VZOUT);
1316 EV_TAB[i][2] = VXOUT;
1317 EV_TAB[i][3] = VYOUT;
1318 EV_TAB[i][4] = 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;
1488 if(NbLam0>0)varntp->
kfis = 20;
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){
1714 if(NbLam0>0)varntp->
kfis = 20;
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;
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);
3487 if (fiss->
ifis > 0) {
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.);
3568 bshell = ecld->
ecgnz[in][iz]- ecld->
vgsld[in][iz];
3569 defbet = ecld->
beta2[in][iz];
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));
3636 bshell = ecld->
ecgnz[ind][izd] - ecld->
vgsld[ind][izd];
3637 defbet = ecld->
beta2[ind][izd];
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));
3689 bshell = ecld->
ecgnz[ind][izd] - ecld->
vgsld[ind][izd];
3690 defbet =ecld->
beta2[ind][izd];
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));
3744 bshell = ecld->
ecgnz[ind][izd] - ecld->
vgsld[ind][izd];
3745 defbet = ecld->
beta2[ind][izd];
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));
3800 bshell = ecld->
ecgnz[ind][izd] - ecld->
vgsld[ind][izd];
3801 defbet = ecld->
beta2[ind][izd];
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));
3856 bshell = ecld->
ecgnz[ind][izd] - ecld->
vgsld[ind][izd];
3857 defbet = ecld->
beta2[ind][izd];
3859 iinert = 0.4 * 931.49 * 1.16*1.16 * std::pow(a-4.,5.0/3.0)*(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));
3912 bshell = ecld->
ecgnz[ind][izd] - ecld->
vgsld[ind][izd];
3913 defbet = ecld->
beta2[ind][izd];
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));
3970 bshell = ecld->
ecgnz[ind][izd] - ecld->
vgsld[ind][izd];
3971 defbet = ecld->
beta2[ind][izd];
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.);
4226 gfactor = 1.+gamma*ecld->
ecgnz[in][iz];
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.);
4238 if(gammaemission==1){
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);
4322 cf =
cram(fiss->
bet,homegasp);
4325 fomega_gs(a,zprf,&k1,&omegags,&homegags);
4326 tauc=
tau(fiss->
bet,homegags,ef,ft);
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;
5457 for(
G4int init_i = 0; init_i < 7; init_i++) {
5461 for(
G4int init_i = 0; init_i < 10; init_i++) {
5465 G4double a = 0.0, z = 0.0, amin = 0.0, amax = 0.0, amin2 = 0.0;
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];
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;
8165 G4double beta=0.,beta1=0.,beta2=0.;
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 ) -
8626 2.0 *
frldm( Zsymm, Nsymm, beta ) -
8627 ecoul( Zsymm, Nsymm, beta, Zsymm, Nsymm, beta, 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;
8910 e1 =
gausshaz(0,E1exc,E1exc_sigma);
8911 if(e1<0.)
goto fis987;
8914 e2 =
gausshaz(0,E2exc,E2exc_sigma);
8915 if(e2<0.)
goto fis988;
8917 (*NbLam0_par) = NbLam0;
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;