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++)
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;
2189#ifdef ABLAXX_IN_GEANT4_MODE
2194 if(dataInterface->
readData() ==
true) {
2195 if(verboseLevel > 0) {
2203 for(
G4int z = 0; z < 99; z++) {
2204 for(
G4int n = 0; n < 154; n++) {
2205 ecld->
ecfnz[n][z] = 0.e0;
2210 ecld->
rms[n][z] = dataInterface->
getRms(n,z);
2214 for(
G4int z = 0; z < 137; z++){
2215 for(
G4int n = 0; n < 251; n++){
2221 for(
G4int z = 0; z < 500; z++) {
2222 for(
G4int a = 0; a < 500; a++) {
2223 pace->
dm[z][a] = dataInterface->
getPace2(z,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++){
2250 masses->
bind[j][i]=0.;
2254 ec2sub->
ecnz[j][i] = 0.0;
2256 masses->
bind[j][i] = dataInterface->
getMexp(j,i)-MP*i -MN*j;
2257 ecld->
vgsld[j][i]=0.;
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));
2287 ec2sub->
ecnz[j][i] = dataInterface->
getMexp(j,i) - (mfrldm[j][i] - e0);
2289 ecld->
vgsld[j][i] =
max(0.0,ec2sub->
ecnz[j][i] - delta_tot);
2297 delete dataInterface;
2311 if(fiss->
zt<83 && fiss->
zt>56){
2346 if(fiss->
zt<84 && fiss->
zt>56)
2383 T_freeze_out_in = -6.5;
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);
2451 (*el) = (*el) + ec2sub->
ecnz[a1-z1][z1];
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);
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;
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;
4574 BSHELLCT = ecld->
ecgnz[in][iz];
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;
5019 G4double pi = 3.141592653589793238e0;
5068 if(flag==1){
goto eflmac311;}
5071 if(masses->
mexpiop[in][iz]==1){
5072 return masses->
bind[in][iz];
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) {
5280 if(verboseLevel > 2) {
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;
5370 RMAX = 1.1 * (ecld->
rms[A1-Z1][Z1]+ecld->
rms[A2-Z2][Z2]) + 2.8;
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++) {
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];
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;
5777 for(
G4int n=0;n<itmax;n++){
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) {
5910 fpace2=pace->
dm[ii][jj];
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;
6004 varntp->
ntrack = IMULTBU + IEV_TAB;
6008 for(
G4int i=0;i<IMULTBU;i++){
6018 varntp->
zvv[intp] = iz;
6019 varntp->
avv[intp] = ia;
6020 varntp->
svv[intp] = -1*is;
6023 G4double v2 = BU_TAB[i][4]*BU_TAB[i][4]+BU_TAB[i][5]*BU_TAB[i][5]+BU_TAB[i][6]*BU_TAB[i][6];
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);
6027 varntp->
pxlab[intp] = etot * BU_TAB[i][4] / c;
6028 varntp->
pylab[intp] = etot * BU_TAB[i][5] / c;
6029 varntp->
pzlab[intp] = etot * BU_TAB[i][6] / c;
6030 varntp->
enerj[intp] = etot - avvmass;
6035 for(
G4int i=0;i<IEV_TAB;i++){
6039 G4int is = EV_TAB[i][5];
6044 varntp->
zvv[intp] = iz;
6045 varntp->
avv[intp] = ia;
6046 varntp->
svv[intp] = -1*is;
6050 G4double v2 = EV_TAB[i][2]*EV_TAB[i][2]+EV_TAB[i][3]*EV_TAB[i][3]+EV_TAB[i][4]*EV_TAB[i][4];
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);
6054 varntp->
pxlab[intp] = etot * EV_TAB[i][2] / c;
6055 varntp->
pylab[intp] = etot * EV_TAB[i][3] / c;
6056 varntp->
pzlab[intp] = etot * EV_TAB[i][4] / c;
6057 varntp->
enerj[intp] = etot - avvmass;
6059 varntp->
zvv[intp] = 0;
6060 varntp->
avv[intp] = 1;
6061 varntp->
svv[intp] = -1;
6064 G4double v2 = EV_TAB[i][2]*EV_TAB[i][2]+EV_TAB[i][3]*EV_TAB[i][3]+EV_TAB[i][4]*EV_TAB[i][4];
6065 G4double gamma = std::sqrt(1.0 - v2 / (c*c));
6068 varntp->
pxlab[intp] = etot * EV_TAB[i][2] / c;
6069 varntp->
pylab[intp] = etot * EV_TAB[i][3] / c;
6070 varntp->
pzlab[intp] = etot * EV_TAB[i][4] / c;
6071 varntp->
enerj[intp] = etot - avvmass;
6073 varntp->
zvv[intp] = iz;
6074 varntp->
avv[intp] = ia;
6075 varntp->
svv[intp] = 0;
6079 varntp->
pxlab[intp] = EV_TAB[i][2];
6080 varntp->
pylab[intp] = EV_TAB[i][3];
6081 varntp->
pzlab[intp] = EV_TAB[i][4];
6082 varntp->
enerj[intp] = std::sqrt(EV_TAB[i][2]*EV_TAB[i][2]+EV_TAB[i][3]*EV_TAB[i][3]+EV_TAB[i][4]*EV_TAB[i][4]);
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;
6464 MU = (
A - ap) * ap /
A;
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));
6488 G4double ALPHA2 = std::sqrt(5.0/(4.0*PI))*ecld->
beta2[IN][IZ];
6489 G4double ALPHA4 = std::sqrt(9.0/(4.0*PI))*ecld->
beta4[IN][IZ];
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;
6543 DEFO_INIT = std::sqrt(5.0/(4.0*PI))*ecld->
beta2[fiss->
at-fiss->
zt][fiss->
zt];
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;
7076 BSHELL = ecld->
ecgnz[in][iz]- ecld->
vgsld[in][iz];
7077 DEFBET = ecld->
beta2[in][iz];
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;
7711 EV_TAB_SSC[IEV_TAB_SSC][0] = 0.;
7712 EV_TAB_SSC[IEV_TAB_SSC][1] = -2.;
7713 EV_TAB_SSC[IEV_TAB_SSC][5] = 1.;
7715 EV_TAB_SSC[IEV_TAB_SSC][0] = zmoins;
7716 EV_TAB_SSC[IEV_TAB_SSC][1] = amoins;
7717 EV_TAB_SSC[IEV_TAB_SSC][5] = 0.;
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);
7733 EV_TAB_SSC[IEV_TAB_SSC][2] = c * pc * xcv / ETOT_LP;
7734 EV_TAB_SSC[IEV_TAB_SSC][3] = c * pc * ycv / ETOT_LP;
7735 EV_TAB_SSC[IEV_TAB_SSC][4] = c * pc * zcv / ETOT_LP;
7738 EV_TAB_SSC[IEV_TAB_SSC][2] = pc * xcv;
7739 EV_TAB_SSC[IEV_TAB_SSC][3] = pc * ycv;
7740 EV_TAB_SSC[IEV_TAB_SSC][4] = pc * zcv;
7743 EV_TAB_SSC[IEV_TAB_SSC][2],EV_TAB_SSC[IEV_TAB_SSC][3],
7744 EV_TAB_SSC[IEV_TAB_SSC][4],
7745 &VXOUT,&VYOUT,&VZOUT);
7746 EV_TAB_SSC[IEV_TAB_SSC][2] = VXOUT;
7747 EV_TAB_SSC[IEV_TAB_SSC][3] = VYOUT;
7748 EV_TAB_SSC[IEV_TAB_SSC][4] = VZOUT;
7752 G4double v2 = std::pow(EV_TAB_SSC[IEV_TAB_SSC][2],2.) +
7753 std::pow(EV_TAB_SSC[IEV_TAB_SSC][3],2.) +
7754 std::pow(EV_TAB_SSC[IEV_TAB_SSC][4],2.);
7755 G4double gamma = 1.0/std::sqrt(1.0 - v2 / (c*c));
7756 G4double etot_lp = amoins*mu * gamma;
7757 pxeva = pxeva - EV_TAB_SSC[IEV_TAB_SSC][2] * etot_lp / c;
7758 pyeva = pyeva - EV_TAB_SSC[IEV_TAB_SSC][3] * etot_lp / c;
7759 pleva = pleva - EV_TAB_SSC[IEV_TAB_SSC][4] * etot_lp / c;
7762 pxeva = pxeva - EV_TAB_SSC[IEV_TAB_SSC][2];
7763 pyeva = pyeva - EV_TAB_SSC[IEV_TAB_SSC][3];
7764 pleva = pleva - EV_TAB_SSC[IEV_TAB_SSC][4];
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;
7773 IEV_TAB_SSC = IEV_TAB_SSC +1;
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)
return 2.04;
7798 else if(Z==2 &&
A==4)
return 2.39;
7799 else if(Z==2 &&
A==5)
return 3.12;
7800 else if(Z==2 &&
A==6)
return 4.18;
7801 else if(Z==2 &&
A==7)
return 5.23;
7802 else if(Z==2 &&
A==8)
return 7.16;
7803 else if(Z==3 &&
A==6)
return 4.50;
7804 else if(Z==3 &&
A==7)
return 5.58;
7805 else if(Z==3 &&
A==8)
return 6.80;
7806 else if(Z==3 &&
A==9)
return 8.50;
7807 else if(Z==4 &&
A==7)
return 5.16;
7808 else if(Z==4 &&
A==8)
return 6.84;
7809 else if(Z==4 &&
A==9)
return 6.71;
7810 else if(Z==4 &&
A==10)
return 9.11;
7811 else if(Z==5 &&
A==9)
return 8.29;
7812 else if(Z==5 &&
A==10)
return 8.89;
7813 else if(Z==5 &&
A==11)
return 10.24;
7814 else if(Z==5 &&
A==12)
return 11.37;
7815 else if(Z==6 &&
A==12)
return 10.76;
7816 else if(Z==6 &&
A==13)
return 11.69;
7817 else if(Z==6 &&
A==14)
return 12.17;
7818 else if(Z==14 &&
A==28)
return 16.0;
7819 else if(Z==39 &&
A==89)
return 22.1;
7820 else if(Z==57 &&
A==139)
return 23.8;
7821 else if(Z==82 &&
A==208)
return 26.5;
7832 if(
A<2 || Z<2)
return 0.;
7842 if(
mod(N,2) == 1 &&
mod(Z,2) == 1)
D = -12./std::sqrt(
A);
7843 if(
mod(N,2) == 0 &&
mod(Z,2) == 0)
D = 12./std::sqrt(
A);
7845 G4double deltanew = (1.-std::exp(-1.*
A/c))*
D;
7847 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));
7851void 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)
7860 e =
dmin1(SBHE,SN,e);
7861 e =
dmin1(SBHE,SBA,e);
7882 *ECP = (-1.0)*SP + BP;
7899 *ECD = (-1.0)*SD + BD;
7916 *ECT = (-1.0)*ST + BT;
7933 *ECHE= (-1.0)*SHE + BHE;
7951 *ECA = (-1.0)*SA + BA;
8033 Delta_U1_shell_max = -2.45;
8039 Delta_U2_shell = -2.45;
8049 G4double Fwidth_asymm1,Fwidth_asymm2,Fwidth_symm;
8051 Fwidth_asymm1 = 0.65;
8052 Fwidth_asymm2 = 0.65;
8094 Friction_factor = 1.0;
8097 G4double cN_asymm1_shell, cN_asymm2_shell;
8098 G4double gamma,gamma_heavy1,gamma_heavy2;
8114 G4double Sasymm1=0.,Sasymm2=0.,Ssymm=0.,Ysum=0.,Yasymm=0.;
8116 G4double wNasymm1_saddle, wNasymm2_saddle, wNsymm_saddle;
8117 G4double wNasymm2_scission, wNsymm_scission;
8118 G4double wNasymm1, wNasymm2, wNsymm;
8141 G4double A_levdens_heavy1,A_levdens_heavy2;
8145 G4double epsilon_1_saddle,epsilon0_1_saddle;
8146 G4double epsilon_2_saddle,epsilon0_2_saddle,epsilon_symm_saddle;
8151 G4double E_eff1_saddle,E_eff2_saddle;
8152 G4double Epot0_mode1_saddle,Epot0_mode2_saddle,Epot0_symm_saddle;
8153 G4double Epot_mode1_saddle,Epot_mode2_saddle,Epot_symm_saddle;
8154 G4double E_defo,E_defo1,E_defo2,E_scission_pre,E_scission_post;
8160 G4double MassCurv_scis, MassCurv_sadd;
8162 G4double Nheavy1_shell,Nheavy2_shell;
8167 G4double Z_scission,N_scission,A_scission;
8171 G4double DSN132,Delta_U1_shell,E_eff0_saddle;
8172 G4int NbLam0= (*NbLam0_par);
8177 cN_asymm1_shell = 0.700 * N/Z;
8178 cN_asymm2_shell = 0.040 * N/Z;
8182 DSN132 = Nheavy1_in - N/Z * Zheavy1_in;
8183 Aheavy1 = Nheavy1_in + Zheavy1_in + 0.340 * DSN132;
8192 Delta_U1_shell = Delta_U1_shell_max + U1NZ_SLOPE * std::abs(DSN132);
8193 Delta_U1_shell =
min(0.,Delta_U1_shell);
8197 Nheavy1 = N/
A * Aheavy1;
8198 Aheavy2 = Nheavy2 *
A/N;
8202 A_levdens =
A / xLevdens;
8203 gamma = A_levdens / (0.40 * std::pow(
A,1.3333)) * FGAMMA;
8204 A_levdens_heavy1 = Aheavy1 / xLevdens;
8205 gamma_heavy1 = A_levdens_heavy1 / (0.40 * std::pow(Aheavy1,1.3333)) * FGAMMA * FGAMMA1;
8206 A_levdens_heavy2 = Aheavy2 / xLevdens;
8207 gamma_heavy2 = A_levdens_heavy2 / (0.40 * std::pow(Aheavy2,1.3333)) * FGAMMA;
8211 E_saddle_scission = (-24. + 0.02227 * Z*Z/std::pow(
A,0.33333))*Friction_factor;
8212 E_saddle_scission =
max( 0.0, E_saddle_scission );
8218 Z2_over_A_eff = Z*Z/
A;
8220 if( Z2_over_A_eff< 34.0 )
8221 MassCurv_scis = std::pow(10., -1.093364 + 0.082933 * Z2_over_A_eff - 0.0002602 * Z2_over_A_eff*Z2_over_A_eff);
8223 MassCurv_scis = std::pow(10., 3.053536 - 0.056477 * Z2_over_A_eff+ 0.0002454 * Z2_over_A_eff*Z2_over_A_eff );
8228 MassCurv_sadd = X_s2s * MassCurv_scis;
8230 cN_symm = 8.0 / std::pow(N,2.) * MassCurv_scis;
8231 cN_symm_sadd = 8.0 / std::pow(N,2.) * MassCurv_sadd;
8233 Nheavy1_shell = Nheavy1;
8236 Nheavy1_eff = (cN_symm_sadd*Nsymm + cN_asymm1_shell *
8237 Uwash(E/
A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1) *
8241 Uwash(E/
A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1));
8243 Nheavy1_eff = (cN_symm_sadd*Nsymm +
8244 cN_asymm1_shell*Nheavy1_shell)
8249 Nheavy2_NZ = Nheavy2;
8250 Nheavy2_shell = Nheavy2_NZ;
8252 Nheavy2_eff = (cN_symm_sadd*Nsymm +
8254 Uwash(E/
A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2) *
8258 Uwash(E/
A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2));
8260 Nheavy2_eff = (cN_symm_sadd*Nsymm +
8261 cN_asymm2_shell*Nheavy2_shell)
8265 Delta_U1 = Delta_U1_shell + (Nheavy1_shell - Nheavy1_eff)*(Nheavy1_shell - Nheavy1_eff) * cN_asymm1_shell;
8266 Delta_U1 =
min(Delta_U1,0.0);
8267 Delta_U2 = Delta_U2_shell + (Nheavy2_shell - Nheavy2_eff)*(Nheavy2_shell - Nheavy2_eff) * cN_asymm2_shell;
8268 Delta_U2 =
min(Delta_U2,0.0);
8272 Epot0_mode1_saddle = (Nheavy1_eff-Nsymm)*(Nheavy1_eff-Nsymm) * cN_symm_sadd;
8273 Epot0_mode2_saddle = (Nheavy2_eff-Nsymm)*(Nheavy2_eff-Nsymm) * cN_symm_sadd;
8274 Epot0_symm_saddle = 0.0;
8278 Epot_mode1_saddle = Epot0_mode1_saddle + Delta_U1;
8279 Epot_mode2_saddle = Epot0_mode2_saddle + Delta_U2;
8280 Epot_symm_saddle = Epot0_symm_saddle;
8283 dUeff =
min( Epot_mode1_saddle, Epot_mode2_saddle);
8284 dUeff =
min( dUeff, Epot_symm_saddle);
8285 dUeff = dUeff - Epot_symm_saddle;
8294 epsilon0_1_saddle = Eld - Epot0_mode1_saddle;
8295 epsilon0_2_saddle = Eld - Epot0_mode2_saddle;
8298 epsilon_1_saddle = Eld - Epot_mode1_saddle;
8299 epsilon_2_saddle = Eld - Epot_mode2_saddle;
8301 epsilon_symm_saddle = Eld - Epot_symm_saddle;
8304 eexc1_saddle = epsilon_1_saddle;
8305 eexc2_saddle = epsilon_2_saddle;
8308 EEXC_MAX =
max( eexc1_saddle, eexc2_saddle);
8309 EEXC_MAX =
max( EEXC_MAX, Eld);
8312 epsilon_1_scission = Eld + E_saddle_scission - Epot_mode1_saddle;
8313 epsilon_2_scission = Eld + E_saddle_scission - Epot_mode2_saddle;
8316 epsilon_symm_scission = Eld + E_saddle_scission - Epot_symm_saddle;
8319 E_eff1_saddle = epsilon0_1_saddle - Delta_U1 *
8320 Uwash(epsilon_1_saddle/
A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1);
8322 if( E_eff1_saddle < A_levdens * hbom1*hbom1)
8323 E_eff1_saddle = A_levdens * hbom1*hbom1;
8326 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_eff1_saddle) /
8328 Uwash(epsilon_1_saddle/
A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)+
8331 E_eff2_saddle = epsilon0_2_saddle -
8333 Uwash(epsilon_2_saddle/
A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2);
8335 if(E_eff2_saddle < A_levdens * hbom2*hbom2)
8336 E_eff2_saddle = A_levdens * hbom2*hbom2;
8339 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_eff2_saddle) /
8341 Uwash(epsilon_2_saddle/
A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)+
8344 E_eff0_saddle = epsilon_symm_saddle;
8345 if(E_eff0_saddle < A_levdens * hbom3*hbom3)
8346 E_eff0_saddle = A_levdens * hbom3*hbom3;
8349 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_eff0_saddle) /
8352 if(epsilon_symm_scission > 0.0 ){
8353 E_HELP =
max(E_saddle_scission,epsilon_symm_scission);
8355 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*(E_HELP)) /
8359 std::sqrt(0.50 * std::sqrt(1.0/A_levdens*E_saddle_scission) /
8366 if( E_saddle_scission == 0.0 ){
8367 wNasymm1_scission = wNasymm1_saddle;
8368 wNasymm2_scission = wNasymm2_saddle;
8370 if( Nheavy1_eff > 75.0 ){
8371 wNasymm1_scission = std::sqrt(21.0)*N/
A;
8372 wNasymm2_scission =
max( 12.8 - 1.0 *(92.0 - Nheavy2_eff),1.0)*N/
A;
8375 wNasymm1_scission = wNasymm1_saddle;
8376 wNasymm2_scission = wNasymm2_saddle;
8380 wNasymm1_scission =
max( wNasymm1_scission, wNasymm1_saddle );
8381 wNasymm2_scission =
max( wNasymm2_scission, wNasymm2_saddle );
8383 wNasymm1 = wNasymm1_scission * Fwidth_asymm1;
8384 wNasymm2 = wNasymm2_scission * Fwidth_asymm2;
8385 wNsymm = wNsymm_scission * Fwidth_symm;
8388 Aheavy1_eff = Nheavy1_eff *
A/N;
8389 Aheavy2_eff = Nheavy2_eff *
A/N;
8391 A_levdens_heavy1 = Aheavy1_eff / xLevdens;
8392 A_levdens_heavy2 = Aheavy2_eff / xLevdens;
8393 gamma_heavy1 = A_levdens_heavy1 / (0.40 * std::pow(Aheavy1_eff,1.3333)) * FGAMMA * FGAMMA1;
8394 gamma_heavy2 = A_levdens_heavy2 / (0.40 * std::pow(Aheavy2_eff,1.3333)) * FGAMMA;
8396 if( epsilon_symm_saddle < A_levdens * hbom3*hbom3)
8397 Ssymm = 2.0 * std::sqrt(A_levdens*A_levdens * hbom3*hbom3) +
8398 (epsilon_symm_saddle - A_levdens * hbom3*hbom3)/hbom3;
8400 Ssymm = 2.0 * std::sqrt(A_levdens*epsilon_symm_saddle);
8404 if( epsilon0_1_saddle < A_levdens * hbom1*hbom1 )
8405 Ssymm_mode1 = 2.0 * std::sqrt(A_levdens*A_levdens * hbom1*hbom1) +
8406 (epsilon0_1_saddle - A_levdens * hbom1*hbom1)/hbom1;
8408 Ssymm_mode1 = 2.0 * std::sqrt( A_levdens*epsilon0_1_saddle );
8410 if( epsilon0_2_saddle < A_levdens * hbom2*hbom2 )
8411 Ssymm_mode2 = 2.0 * std::sqrt(A_levdens*A_levdens * hbom2*hbom2) +
8412 (epsilon0_2_saddle - A_levdens * hbom2*hbom2)/hbom2;
8414 Ssymm_mode2 = 2.0 * std::sqrt(A_levdens*epsilon0_2_saddle);
8417 if( epsilon0_1_saddle -
8419 Uwash(epsilon_1_saddle/
A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)
8420 < A_levdens * hbom1*hbom1 )
8421 Sasymm1 = 2.0 * std::sqrt( A_levdens*A_levdens * hbom1*hbom1 ) +
8422 (epsilon0_1_saddle - Delta_U1 *
8423 Uwash(epsilon_1_saddle/
A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)
8424 - A_levdens * hbom1*hbom1)/hbom1;
8426 Sasymm1 = 2.0 *std::sqrt( A_levdens*(epsilon0_1_saddle - Delta_U1 *
8427 Uwash(epsilon_1_saddle/
A*Aheavy1,Ecrit,FREDSHELL,gamma_heavy1)));
8429 if( epsilon0_2_saddle -
8431 Uwash(epsilon_2_saddle/
A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)
8432 < A_levdens * hbom2*hbom2 )
8433 Sasymm2 = 2.0 * std::sqrt( A_levdens*A_levdens * hbom2*hbom2 ) +
8434 (epsilon0_1_saddle-Delta_U1 *
8435 Uwash(epsilon_2_saddle/
A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)
8436 - A_levdens * hbom2*hbom2)/hbom2;
8439 std::sqrt( A_levdens*(epsilon0_2_saddle - Delta_U2 *
8440 Uwash(epsilon_2_saddle/
A*Aheavy2,Ecrit,FREDSHELL,gamma_heavy2)));
8442 Yasymm1 = ( std::exp(Sasymm1 - Ssymm) - std::exp(Ssymm_mode1 - Ssymm) ) *
8443 wNasymm1_saddle / wNsymm_saddle * 2.0;
8445 Yasymm2 = ( std::exp(Sasymm2 - Ssymm) - std::exp(Ssymm_mode2 - Ssymm) ) *
8446 wNasymm2_saddle / wNsymm_saddle * 2.0;
8448 Ysum = Ysymm + Yasymm1 + Yasymm2;
8451 Ysymm = Ysymm / Ysum;
8452 Yasymm1 = Yasymm1 / Ysum;
8453 Yasymm2 = Yasymm2 / Ysum;
8454 Yasymm = Yasymm1 + Yasymm2;
8460 if( (epsilon_symm_saddle < epsilon_1_saddle) &&
8461 (epsilon_symm_saddle < epsilon_2_saddle) )
8464 if( epsilon_1_saddle < epsilon_2_saddle )
8472 r_e_o = std::pow(10.0,-0.0170 * (E_saddle_scission + Eld)*(E_saddle_scission + Eld));
8490 if( rmode < Yasymm1 )
8493 if( (rmode > Yasymm1) && (rmode < Yasymm) )
8502 N1mean = Nheavy1_eff;
8506 N1mean = Nheavy2_eff;
8522 while( N1r < 5.0 || N2r < 5.0 ){
8541 E_scission_pre =
max( epsilon_1_scission, 1.0 );
8543 if( N1mean > N*0.50 ){
8553 E_scission_pre =
max( epsilon_2_scission, 1.0 );
8554 if( N1mean > N*0.50 ){
8555 beta1 = (N1r - 92.0) * 0.030 + 0.60;
8560 beta1 =
max(beta1,beta1gs);
8561 beta2 = 1.0 - beta1;
8562 beta2 =
max(beta2,beta2gs);
8568 beta2 = (N2r -92.0) * 0.030 + 0.60;
8569 beta2 =
max(beta2,beta2gs);
8570 beta1 = 1.0 - beta2;
8571 beta1 =
max(beta1,beta1gs);
8586 beta =
max(0.177963+0.0153241*Zsymm-1.62037e-4*Zsymm*Zsymm,betags);
8587 beta1 =
max(0.177963+0.0153241*Z1UCD-1.62037e-4*Z1UCD*Z1UCD,beta1gs);
8588 beta2 =
max(0.177963+0.0153241*Z2UCD-1.62037e-4*Z2UCD*Z2UCD,beta2gs);
8590 E_asym =
frldm( Z1UCD, N1r, beta1 ) +
8591 frldm( Z2UCD, N2r, beta2 ) +
8592 ecoul( Z1UCD, N1r, beta1, Z2UCD, N2r, beta2, 2.0 ) -
8593 2.0 *
frldm( Zsymm, Nsymm, beta ) -
8594 ecoul( Zsymm, Nsymm, beta, Zsymm, Nsymm, beta, 2.0 );
8595 E_scission_pre =
max( epsilon_symm_scission - E_asym, 1. );
8604 if(E_scission_pre>5. && NbLam0<1){
8606 &A_scission,&Z_scission,vx_eva_sc,vy_eva_sc,vz_eva_sc,&NbLam0);
8607 N_scission = A_scission - Z_scission;
8611 E_scission_post = E_scission_pre;
8612 N_scission = A_scission - Z_scission;
8618 N1r = N1r * N_scission / N;
8619 N2r = N2r * N_scission / N;
8620 Z1UCD = Z1UCD * Z_scission / Z;
8621 Z2UCD = Z2UCD * Z_scission / Z;
8659 CZ = (
frldm( Z1UCD-1.0, N1r+1.0, beta1 ) +
8660 frldm( Z2UCD+1.0, N2r-1.0, beta2 ) +
8661 frldm( Z1UCD+1.0, N1r-1.0, beta1 ) +
8662 frldm( Z2UCD-1.0, N2r+1.0, beta2 ) +
8663 ecoul( Z1UCD-1.0, N1r+1.0, beta1,
8664 Z2UCD+1.0, N2r-1.0, beta2, 2.0) +
8665 ecoul( Z1UCD+1.0, N1r-1.0, beta1,
8666 Z2UCD-1.0, N2r+1.0, beta2, 2.0) -
8667 2.0*
ecoul( Z1UCD, N1r, beta1, Z2UCD, N2r, beta2, 2.0) -
8668 2.0*
frldm( Z1UCD, N1r, beta1 ) -
8669 2.0*
frldm( Z2UCD, N2r, beta2) ) * 0.50;
8671 if(1.0/A_levdens*E_scission_post < 0.0)
8672 std::cout <<
"DSQRT 1 < 0" << A_levdens <<
" " << E_scission_post << std::endl;
8674 if(0.50 * std::sqrt(1.0/A_levdens*E_scission_post) / CZ < 0.0){
8675 std::cout <<
"DSQRT 2 < 0 " << CZ << std::endl;
8676 std::cout <<
"This event was not considered" << std::endl;
8680 ZA1width = std::sqrt(0.5*std::sqrt(1.0/A_levdens*E_scission_post)/CZ);
8691 ZA1width =
max(ZA1width,sigZmin);
8693 if(imode == 1 && cpol1 != 0.0){
8697 Z1rr = Z1UCD - cpol1 * A_scission/N_scission;
8703 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED" << std::endl;
8706 if ((
utilabs(Z1rr - Z1r) > 3.0*ZA1width) || Z1r<1.0)
goto fiss2801;
8709 if( imode == 2 && cpol2 != 0.0 ){
8713 Z1rr = Z1UCD - cpol2 * A_scission/N_scission;
8718 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED" << std::endl;
8721 if( (
utilabs(Z1rr - Z1r) > 3.0*ZA1width) || Z1r < 1.0 )
goto fiss2802;
8729 re1 =
frldm( Z1UCD-1.0, N1r+1.0, beta1 ) +
8730 frldm( Z2UCD+1.0, N2r-1.0, beta2 ) +
8731 ecoul( Z1UCD-1.0, N1r+1.0, beta1,
8732 Z2UCD+1.0, N2r-1.0, beta2, d );
8733 re2 =
frldm( Z1UCD, N1r, beta1) +
8734 frldm( Z2UCD, N2r, beta2 ) +
8735 ecoul( Z1UCD, N1r, beta1,
8736 Z2UCD, N2r, beta2, d );
8737 re3 =
frldm( Z1UCD+1.0, N1r-1.0, beta1 ) +
8738 frldm( Z2UCD-1.0, N2r+1.0, beta2 ) +
8739 ecoul( Z1UCD+1.0, N1r-1.0, beta1,
8740 Z2UCD-1.0, N2r+1.0, beta2, d );
8741 eps2 = ( re1 - 2.0*re2 + re3 ) / 2.0;
8742 eps1 = ( re3 - re1 ) / 2.0;
8743 DN1_POL = -eps1 / ( 2.0 * eps2 );
8745 Z1rr = Z1UCD + DN1_POL;
8750 DN1_POL = DN1_POL - 0.6 *
Uwash(E_scission_post,Ecrit,FREDSHELL,gamma);
8751 Z1rr = Z1UCD + DN1_POL;
8752 if ( Z1rr < 50. ) Z1rr = 50.0;
8754 DN1_POL = DN1_POL + 0.60 *
Uwash(E_scission_post,Ecrit,FREDSHELL,gamma);
8755 Z1rr = Z1UCD + DN1_POL;
8756 if ( Z1rr > 50.0 ) Z1rr = 50.0;
8766 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING Z1R IN PROFI.FOR. A VALUE WILL BE FORCED" << std::endl;
8770 if( (
utilabs(Z1rr - Z1r) > 3.0*ZA1width) || (Z1r < 1.0) )
goto fiss2803;
8782 z2 =
dint( Z_scission ) - z1;
8784 N2 =
dint( N_scission ) - N1;
8788 if( (z1 < 0) || (z2 < 0) || (a1 < 0) || (a2 < 0) ){
8789 std::cout <<
" -------------------------------" << std::endl;
8790 std::cout <<
" Z, A, N : " << Z <<
" " <<
A <<
" " << N << std::endl;
8791 std::cout << z1 <<
" " << z2 <<
" " << a1 <<
" " << a2 << std::endl;
8792 std::cout << E_scission_post <<
" " << A_levdens <<
" " << CZ << std::endl;
8794 std::cout <<
" -------------------------------" << std::endl;
8803 if( N1mean > N*0.50 ){
8807 if(beta2< beta2gs) beta2 = beta2gs;
8808 E1exc = E_scission_pre * a1 /
A + E_defo;
8809 E_defo =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8810 E2exc = E_scission_pre * a2 /
A + E_defo;
8814 if(beta1< beta1gs) beta1 = beta1gs;
8815 E_defo =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8816 E1exc = E_scission_pre * a1 /
A + E_defo;
8818 E2exc = E_scission_pre * a2 /
A + E_defo;
8825 if( N1mean > N*0.5 ){
8828 if(beta1< beta1gs) beta1 = beta1gs;
8829 E_defo =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8830 E1exc = E_scission_pre * a1 /
A + E_defo;
8832 if(beta2< beta2gs) beta2 = beta2gs;
8833 E_defo =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8834 E2exc = E_scission_pre * a2 /
A + E_defo;
8838 if(beta2< beta2gs) beta2 = beta2gs;
8839 E_defo =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8840 E2exc = E_scission_pre * a2 /
A + E_defo;
8842 if(beta1< beta1gs) beta1 = beta1gs;
8843 E_defo =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8844 E1exc = E_scission_pre * a1 /
A + E_defo;
8851 if(beta1< beta1gs) beta1 = beta1gs;
8853 if(beta2< beta2gs) beta2 = beta2gs;
8854 E_defo1 =
frldm( z1, N1, beta1 ) -
frldm( z1, N1, beta1gs );
8855 E_defo2 =
frldm( z2, N2, beta2 ) -
frldm( z2, N2, beta2gs );
8856 E1exc = E_scission_pre * a1 /
A + E_defo1;
8857 E2exc = E_scission_pre * a2 /
A + E_defo2;
8862 TKER = ( z1 * z2 * 1.440 ) /
8863 ( R0 * std::pow(a1,0.333330) * (1.0 + 2.0/3.0 * beta1 ) +
8864 R0 * std::pow(a2,0.333330) * (1.0 + 2.0/3.0 * beta2 ) + 2.0 );
8866 EkinR1 = TKER * a2 /
A;
8867 EkinR2 = TKER * a1 /
A;
8868 v1 = std::sqrt(EkinR1/a1) * 1.3887;
8869 v2 = std::sqrt(EkinR2/a2) * 1.3887;
8877 e1 =
gausshaz(0,E1exc,E1exc_sigma);
8878 if(e1<0.)
goto fis987;
8881 e2 =
gausshaz(0,E2exc,E2exc_sigma);
8882 if(e2<0.)
goto fis988;
8884 (*NbLam0_par) = NbLam0;
8922 G4double r_in = 0.0, r_rest = 0.0, r_help = 0.0;
8928 r_in = r_origin + 0.5;
8930 if (r_even_odd < 0.001) {
8931 i_out = (
G4int)(r_floor);
8934 r_rest = r_in - r_floor;
8935 r_middle = r_floor + 0.5;
8936 n_floor = (
G4int)(r_floor);
8937 if (n_floor%2 == 0) {
8939 r_help = r_middle + (r_rest - 0.5) * (1.0 - r_even_odd);
8943 r_help = r_middle + (r_rest - 0.5) * (1.0 + r_even_odd);
8945 i_out = (
G4int)(r_help);
8961 G4double xcom = 0.0, xvs = 0.0, xe = 0.0;
8965 alpha = ( std::sqrt(5.0/(4.0*pi)) ) * beta;
8967 xcom = 1.0 - 1.7826 * ((a - 2.0*z)/a)*((a - 2.0*z)/a);
8969 xvs = - xcom * ( 15.4941 * a -
8970 17.9439 * std::pow(a,2.0/3.0) * (1.0+0.4*alpha*alpha) );
8972 xe = z*z * (0.7053/(std::pow(a,1.0/3.0)) * (1.0-0.2*alpha*alpha) - 1.1529/a);
8999 dtot = r0 * ( std::pow((z1+n1),1.0/3.0) * (1.0+0.6666667*beta1)
9000 + std::pow((z2+n2),1.0/3.0) * (1.0+0.6666667*beta2) ) + d;
9001 fecoul = z1 * z2 * 1.44 / dtot;
9013 R_wash = std::exp(-E * Freduction * gamma);
9015 R_wash = std::exp(- Ecrit * Freduction * gamma -(E-Ecrit) * gamma);
9073 G4double z = 0.0, n = 0.0, a = 0.0, av = 0.0, as = 0.0;
9074 G4double a0 = 0.0, c1 = 0.0, c4 = 0.0, b1 = 0.0, b3 = 0.0;
9075 G4double ff = 0.0, ca = 0.0, w = 0.0, efl = 0.0;
9076 G4double r0 = 0.0, kf = 0.0, ks = 0.0;
9077 G4double kv = 0.0, rp = 0.0, ay = 0.0, aden = 0.0, x0 = 0.0, y0 = 0.0;
9078 G4double esq = 0.0, ael = 0.0, i = 0.0;
9079 G4double pi = 3.141592653589793238e0;
9129 c1 = 3.0/5.0*esq/r0;
9130 c4 = 5.0/4.0*std::pow((3.0/(2.0*pi)),(2.0/3.0)) * c1;
9131 kf = std::pow((9.0*pi*z/(4.0*a)),(1.0/3.0))/r0;
9133 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));
9137 x0 = r0 * std::pow(a,(1.0/3.0)) / ay;
9138 y0 = r0 * std::pow(a,(1.0/3.0)) / aden;
9140 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);
9142 b3 = 1.0 - 5.0/std::pow(y0,2) * (1.0 - 15.0/(8.0*y0) + 21.0/(8.0 * std::pow(y0,3))
9143 - 3.0/4.0 * (1.0 + 9.0/(2.0*y0) + 7.0/std::pow(y0,2)
9144 + 7.0/(2.0 * std::pow(y0,3))) * std::exp(-2.0*y0));
9148 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
9149 + 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))
9150 + ff*std::pow(z,2)/a -ca*(n-z) - ael * std::pow(z,(2.39e0));
9156 return eflmacResult;
9161void 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){
9163 G4int INMIN,INMAX,NDIF=0,IMEM;
9164 G4int NEVA=0,PEVA=0;
9172 for(
G4int i=0;i<200;i++){
9173 BU_TAB_TEMP[i][0] = 0.0;
9174 BU_TAB_TEMP[i][1] = 0.0;
9175 BU_TAB_TEMP[i][2] = 0.0;
9176 BU_TAB_TEMP[i][3] = 0.0;
9177 BU_TAB_TEMP[i][4] = 0.0;
9184 if(AFP==0 && ZFP==0){
9188 if((AFP==1 && ZFP==0) ||
9189 (AFP==1 && ZFP==1) ||
9190 (AFP==2 && ZFP==1) ||
9191 (AFP==3 && ZFP==1) ||
9192 (AFP==3 && ZFP==2) ||
9193 (AFP==4 && ZFP==2) ||
9194 (AFP==6 && ZFP==2) ||
9203 if ((AFP-ZFP)==0 && ZFP>1){
9204 for(
G4int I = 0;I<=AFP-2;I++){
9206 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9207 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9208 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9209 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9210 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9211 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9212 *ILOOP = *ILOOP + 1;
9229 for(
G4int I = 1;I<=10; I++){
9241 for(
G4int I = 0;I< IMEM;I++){
9247 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9248 BU_TAB_TEMP[I+1+*ILOOP][0] = 1.0;
9249 BU_TAB_TEMP[I+1+*ILOOP][1] = 1.0;
9250 BU_TAB_TEMP[I+1+*ILOOP][2] = VP2X;
9251 BU_TAB_TEMP[I+1+*ILOOP][3] = VP2Y;
9252 BU_TAB_TEMP[I+1+*ILOOP][4] = VP2Z;
9257 *ILOOP = *ILOOP + IMEM;
9262 NEVA = NDIF - INMAX;
9265 for(
G4int I = 0;I<NEVA;I++){
9271 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9273 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9274 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9275 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9276 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9277 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9278 *ILOOP = *ILOOP + 1;
9285 if ((AFP>=2) && (ZFP==0)){
9286 for(
G4int I = 0;I<= AFP-2;I++){
9290 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9292 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9293 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9294 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9295 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9296 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9297 *ILOOP = *ILOOP + 1;
9309 std::cout <<
"WARNING - BU UNSTABLE: AF < ZF" << std::endl;
9314 if ((AFP>=4) && (ZFP==1)){
9316 for(
G4int I = 0; I<AFP-3;I++){
9320 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9322 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9323 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9324 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9325 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9326 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9327 *ILOOP = *ILOOP + 1;
9339 if ((AFP==4) && (ZFP==3)){
9347 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9349 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9350 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9351 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9352 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9353 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9354 *ILOOP = *ILOOP + 1;
9356 if ((AFP==5) && (ZFP==2)){
9364 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9365 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9366 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9367 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9368 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9369 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9370 *ILOOP = *ILOOP + 1;
9373 if ((AFP==5) && (ZFP==3)){
9381 &(*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;
9390 if ((AFP==6) && (ZFP==4)){
9399 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9400 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9401 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9402 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9403 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9404 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9405 *ILOOP = *ILOOP + 1;
9413 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9414 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9415 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9416 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9417 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9418 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9419 *ILOOP = *ILOOP + 1;
9421 if ((AFP==7)&&(ZFP==2)){
9429 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9430 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9431 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9432 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9433 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9434 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9435 *ILOOP = *ILOOP + 1;
9438 if ((AFP==7) && (ZFP==5)){
9440 for(
int I = 0; I<= AFP-5;I++){
9442 double(AFP-I-1),
double(ZFP-I-1),
9444 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9445 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9446 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9447 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9448 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9449 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9450 *ILOOP = *ILOOP + 1;
9461 if ((AFP==8) && (ZFP==4)){
9468 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9469 BU_TAB_TEMP[*ILOOP][0] = 2.0;
9470 BU_TAB_TEMP[*ILOOP][1] = 4.0;
9471 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9472 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9473 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9474 *ILOOP = *ILOOP + 1;
9476 if ((AFP==8) && (ZFP==6)){
9485 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9486 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9487 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9488 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9489 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9490 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9491 *ILOOP = *ILOOP + 1;
9498 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9499 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9500 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9501 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9502 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9503 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9504 *ILOOP = *ILOOP + 1;
9510 if((AFP==9) && (ZFP==2)){
9519 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9520 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9521 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9522 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9523 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9524 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9525 *ILOOP = *ILOOP + 1;
9531 if((AFP==9) && (ZFP==5)){
9539 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9540 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9541 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9542 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9543 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9544 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9545 *ILOOP = *ILOOP + 1;
9552 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9553 BU_TAB_TEMP[*ILOOP][0] = 2.0;
9554 BU_TAB_TEMP[*ILOOP][1] = 4.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==10) && (ZFP==2)){
9573 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9574 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9575 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9576 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9577 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9578 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9579 *ILOOP = *ILOOP + 1;
9587 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9588 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9589 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9590 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9591 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9592 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9593 *ILOOP = *ILOOP + 1;
9598 if ((AFP==10) && (ZFP==3)){
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;
9617 if ((AFP==10) && (ZFP==7)){
9625 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9626 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9627 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9628 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9629 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9630 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9631 *ILOOP = *ILOOP + 1;
9637 if((AFP==11) && (ZFP==7)){
9645 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9646 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9647 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9648 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9649 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9650 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9651 *ILOOP = *ILOOP + 1;
9656 if ((AFP==12) && (ZFP==8)){
9665 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9666 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9667 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9668 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9669 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9670 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9671 *ILOOP = *ILOOP + 1;
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==15) && (ZFP==9)){
9697 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9698 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9699 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9700 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9701 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9702 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9703 *ILOOP = *ILOOP + 1;
9709 if ((AFP==16) && (ZFP==9)){
9717 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9718 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9719 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9720 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9721 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9722 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9723 *ILOOP = *ILOOP + 1;
9729 if ((AFP==16) && (ZFP==10)){
9737 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9738 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9739 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9740 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9741 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9742 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9743 *ILOOP = *ILOOP + 1;
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;
9761 if((AFP==18) && (ZFP==11)){
9769 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9770 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9771 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9772 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9773 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9774 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9775 *ILOOP = *ILOOP + 1;
9780 if((AFP==19) && (ZFP==11)){
9788 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9789 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9790 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9791 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9792 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9793 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9794 *ILOOP = *ILOOP + 1;
9799 if (ZFP>=4 && (AFP-ZFP)==1){
9804 for(
G4int I = 0;I< NEVA;I++){
9808 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9809 BU_TAB_TEMP[*ILOOP][0] = 0.0;
9810 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9811 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9812 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9813 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9814 *ILOOP = *ILOOP + 1;
9819 for(
G4int I = 0;I<PEVA;I++){
9823 &(*VP1X),&(*VP1Y),&(*VP1Z),&VP2X,&VP2Y,&VP2Z);
9824 BU_TAB_TEMP[*ILOOP][0] = 1.0;
9825 BU_TAB_TEMP[*ILOOP][1] = 1.0;
9826 BU_TAB_TEMP[*ILOOP][2] = VP2X;
9827 BU_TAB_TEMP[*ILOOP][3] = VP2Y;
9828 BU_TAB_TEMP[*ILOOP][4] = VP2Z;
9829 *ILOOP = *ILOOP + 1;
9847void 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){
9850 G4double PX1,PX2,PY1,PY2,PZ1,PZ2,PTOT;
9851 G4double RNDT,CTET1,STET1,RNDP,PHI1,ETOT_P1,ETOT_P2;
9853 G4double vxout=0.,vyout=0.,vzout=0.;
9854 G4int iain,izin,ianew,iznew,inin,innew;
9864 innew = ianew - iznew;
9869 mglms(ain,zin,3,&MASS);
9870 mglms(anew,znew,3,&MASS1);
9871 mglms(ain-anew,zin-znew,3,&MASS2);
9872 ekin_tot = MASS-MASS1-MASS2;
9875 ekin_tot = masses->
massexp[inin][izin]-(masses->
massexp[innew][iznew]+masses->
massexp[inin-innew][izin-iznew]);
9876 if(izin>12)std::cout <<
"*** ZIN > 12 ***" << izin << std::endl;
9879 if( ekin_tot<0.00 ){
9890 EKIN_P1 = ekin_tot*(ain-anew)/ ain;
9891 ETOT_P1 = EKIN_P1 + anew * AMU;
9892 PTOT = anew*AMU*std::sqrt((EKIN_P1/(anew*AMU)+1.0)*(EKIN_P1/(anew*AMU)+1.0)-1.0);
9895 CTET1 = 2.0*RNDT-1.0;
9896 STET1 = std::sqrt(1.0-CTET1*CTET1);
9898 PHI1 = RNDP*2.0*3.141592654;
9899 PX1 = PTOT * STET1*std::cos(PHI1);
9900 PY1 = PTOT * STET1*std::sin(PHI1);
9902 *v1x =
C * PX1 / ETOT_P1;
9903 *v1y =
C * PY1 / ETOT_P1;
9904 *v1z =
C * PZ1 / ETOT_P1;
9905 lorentz_boost(vxin,vyin,vzin,*v1x,*v1y,*v1z,&vxout,&vyout,&vzout);
9913 ETOT_P2 = (ekin_tot - EKIN_P1) + (ain-anew) * AMU;
9914 *v2x =
C * PX2 / ETOT_P2;
9915 *v2y =
C * PY2 / ETOT_P2;
9916 *v2z =
C * PZ2 / ETOT_P2;
9917 lorentz_boost(vxin,vyin,vzin,*v2x,*v2y,*v2z,&vxout,&vyout,&vzout);
9935 G4double GAMMA,VR,
C,CC,DENO,VXNOM,VYNOM,VZNOM;
9946 VR = std::sqrt(VXR*VXR + VYR*VYR + VZR*VZR);
9953 GAMMA = 1.0/std::sqrt(1.0 - VR*VR/CC);
9954 DENO = 1.0 - VXR*VXIN/CC - VYR*VYIN/CC - VZR*VZIN/CC;
9957 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;
9959 *VXOUT = VXNOM / (GAMMA * DENO);
9962 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;
9964 *VYOUT = VYNOM / (GAMMA * DENO);
9967 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;
9969 *VZOUT = VZNOM / (GAMMA * DENO);
9981 G4double EFF1=0.,EFF2=0.,VFF1=0.,VFF2=0.,
9982 AF1=0.,ZF1=0.,AF2=0.,ZF2=0.,
9983 AFF1=0.,ZFF1=0.,AFF2=0.,ZFF2=0.,
9984 vz1_eva=0., vx1_eva=0.,vy1_eva=0.,
9985 vz2_eva=0., vx2_eva=0.,vy2_eva=0.,
9986 vx_eva_sc=0.,vy_eva_sc=0.,vz_eva_sc=0.,
9987 VXOUT=0.,VYOUT=0.,VZOUT=0.,
9988 VX2OUT=0.,VY2OUT=0.,VZ2OUT=0.;
9989 G4int IEV_TAB_FIS=0,IEV_TAB_TEMP=0;
9990 G4double EV_TEMP1[200][6], EV_TEMP2[200][6],mtota=0.;
9991 G4int inttype = 0,inum=0;
9994 G4int NbLam0= (*NbLam0_par);
9996 for(
G4int I1=0;I1<200;I1++)
9997 for(
G4int I2=0;I2<6;I2++){
9998 EV_TEMP[I1][I2] = 0.0;
9999 EV_TEMP1[I1][I2] = 0.0;
10000 EV_TEMP2[I1][I2] = 0.0;
10003 G4double et = EE - JPRF * JPRF * 197. * 197./(2.*0.4*931.*std::pow(AF,5.0/3.0)*1.16*1.16);
10005 fissionDistri(AF,ZF,et,AF1,ZF1,EFF1,VFF1,AF2,ZF2,EFF2,VFF2,
10006 vx_eva_sc,vy_eva_sc,vz_eva_sc,&NbLam0);
10011 G4double pbH = (AF1 - ZF1) / (AF1 - ZF1 + AF2 - ZF2);
10012 for(
G4int i=0;i<NbLam0;i++){
10020 for(
G4int IJ = 0; IJ< IEV_TAB_SSC;IJ++){
10021 EV_TEMP[IJ][0] = EV_TAB_SSC[IJ][0];
10022 EV_TEMP[IJ][1] = EV_TAB_SSC[IJ][1];
10023 EV_TEMP[IJ][2] = EV_TAB_SSC[IJ][2];
10024 EV_TEMP[IJ][3] = EV_TAB_SSC[IJ][3];
10025 EV_TEMP[IJ][4] = EV_TAB_SSC[IJ][4];
10026 EV_TEMP[IJ][5] = EV_TAB_SSC[IJ][5];
10028 IEV_TAB_FIS = IEV_TAB_FIS + IEV_TAB_SSC;
10032 G4double VPERP1 = std::sqrt(VFF1*VFF1 - VZ1_FISSION*VZ1_FISSION);
10034 G4double VX1_FISSION = VPERP1 * std::sin(ALPHA1);
10035 G4double VY1_FISSION = VPERP1 * std::cos(ALPHA1);
10036 G4double VX2_FISSION = - VX1_FISSION / VFF1 * VFF2;
10037 G4double VY2_FISSION = - VY1_FISSION / VFF1 * VFF2;
10038 G4double VZ2_FISSION = - VZ1_FISSION / VFF1 * VFF2;
10041 if( (ZF1<=0.0) || (AF1<=0.0) || (AF1<ZF1) ){
10042 std::cout <<
"F1 unphysical: "<<ZF<<
" "<<AF<<
" "<<EE<<
" "<<ZF1<<
" "<<AF1 << std::endl;
10048 G4int FF11=0, FIMF11=0;
10049 G4double ZIMFF1=0., AIMFF1=0.,TKEIMF1=0.,JPRFOUT=0.;
10051 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);
10053 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
10054 EV_TEMP[IJ+IEV_TAB_FIS][0] = EV_TEMP1[IJ][0];
10055 EV_TEMP[IJ+IEV_TAB_FIS][1] = EV_TEMP1[IJ][1];
10062 EV_TEMP1[IJ][2],EV_TEMP1[IJ][3],EV_TEMP1[IJ][4],
10063 &VXOUT,&VYOUT,&VZOUT);
10066 &VX2OUT,&VY2OUT,&VZ2OUT);
10067 EV_TEMP[IJ+IEV_TAB_FIS][2] = VX2OUT;
10068 EV_TEMP[IJ+IEV_TAB_FIS][3] = VY2OUT;
10069 EV_TEMP[IJ+IEV_TAB_FIS][4] = VZ2OUT;
10072 IEV_TAB_FIS = IEV_TAB_FIS + IEV_TAB_TEMP;
10077 if( (ZF2<=0.0) || (AF2<=0.0) || (AF2<ZF2) ){
10078 std::cout <<
"F2 unphysical: "<<ZF<<
" "<<AF<<
" "<<EE<<
" "<<ZF2<<
" "<<AF2 << std::endl;
10084 G4int FF22=0, FIMF22=0;
10085 G4double ZIMFF2=0., AIMFF2=0.,TKEIMF2=0.,JPRFOUT=0.;
10087 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);
10089 for(
G4int IJ = 0; IJ< IEV_TAB_TEMP;IJ++){
10090 EV_TEMP[IJ+IEV_TAB_FIS][0] = EV_TEMP2[IJ][0];
10091 EV_TEMP[IJ+IEV_TAB_FIS][1] = EV_TEMP2[IJ][1];
10098 EV_TEMP2[IJ][2],EV_TEMP2[IJ][3],EV_TEMP2[IJ][4],
10099 &VXOUT,&VYOUT,&VZOUT);
10102 &VX2OUT,&VY2OUT,&VZ2OUT);
10103 EV_TEMP[IJ+IEV_TAB_FIS][2] = VX2OUT;
10104 EV_TEMP[IJ+IEV_TAB_FIS][3] = VY2OUT;
10105 EV_TEMP[IJ+IEV_TAB_FIS][4] = VZ2OUT;
10108 IEV_TAB_FIS = IEV_TAB_FIS + IEV_TAB_TEMP;
10121 VX1_FISSION,VY1_FISSION,VZ1_FISSION,
10122 &VXOUT,&VYOUT,&VZOUT);
10123 VX1_FISSION = VXOUT;
10124 VY1_FISSION = VYOUT;
10125 VZ1_FISSION = VZOUT;
10127 VX2_FISSION,VY2_FISSION,VZ2_FISSION,
10128 &VXOUT,&VYOUT,&VZOUT);
10129 VX2_FISSION = VXOUT;
10130 VY2_FISSION = VYOUT;
10131 VZ2_FISSION = VZOUT;
10136 (*VX1_FISSION_par) = VX1_FISSION;
10137 (*VY1_FISSION_par) = VY1_FISSION;
10138 (*VZ1_FISSION_par) = VZ1_FISSION;
10139 (*VX_EVA_SC_par)=vx_eva_sc;
10140 (*VY_EVA_SC_par)=vy_eva_sc;
10141 (*VZ_EVA_SC_par)=vz_eva_sc;
10145 (*VX2_FISSION_par) = VX2_FISSION;
10146 (*VY2_FISSION_par) = VY2_FISSION;
10147 (*VZ2_FISSION_par) = VZ2_FISSION;
10148 (*IEV_TAB_FIS_par) = IEV_TAB_FIS;
10149 (*NbLam0_par) = NbLam1 + NbLam2;
10150 if(NbLam0>(NbLam1 + NbLam2))varntp->
kfis = 25;
10157 G4double V_over_V0,R0,RALL,RHAZ,R,TKE,Ekin,V,VPERP,ALPHA1;
10169 RALL = R0 * std::pow(V_over_V0,1.0/3.0) * std::pow(AAL,1.0/3.0);
10171 R = std::pow(RHAZ,1.0/3.0) * RALL;
10172 TKE = 1.44 * Z * ZALL * R*R * (1.0 -
A/AAL)*(1.0 -
A/AAL) / std::pow(RALL,3.0);
10174 Ekin = TKE * (AAL -
A) / AAL;
10176 V = std::sqrt(Ekin/
A) * 1.3887;
10178 VPERP = std::sqrt(V*V - (*VZ)*(*VZ));
10180 *VX = VPERP * std::sin(ALPHA1);
10181 *VY = VPERP * std::cos(ALPHA1);
10207 ix =
G4int(y * 100 + 43543000);
10208 if(
mod(ix,2) == 0) {
10222 l_plus = lambda + 1.;
10226 y=std::pow(
G4AblaRandom::flat()*(std::pow(rxmax,l_plus)-std::pow(rxmin,l_plus))+ std::pow(rxmin,l_plus),1.0/l_plus);
10260 G4double Efermi = 5.0 * SIGMA_0 * SIGMA_0 / (2.0 * 931.4940);
10266 SIGMA_0 = SIGMA_0 * std::pow(V0_over_VBU,1.0/3.0);
10271 GOLDHA_BU = SIGMA_0 * std::sqrt((APRF*(AABRA-APRF))/(AABRA-1.0));
10272 GOLDHA = GOLDHA_BU*std::sqrt(1.0 +
10273 5.0 * PI*PI / 12.0 * (T_freeze_out / Efermi)*(T_freeze_out / Efermi));
10279 GOLDHA_BU = std::sqrt(APRF * T_freeze_out * 931.494 *
10280 (AABRA - APRF) / AABRA);
10281 GOLDHA = GOLDHA_BU;
10286 GOLDHA = SIGMA_0 * std::sqrt((APRF*(AABRA-APRF))/(AABRA-1.0));
10294 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING PX IN Rn07.FOR. A VALUE WILL BE FORCED." << std::endl;
10295 *PX = (AABRA-1.0)*931.4940;
10297 if(std::abs(*PX)>= AABRA*931.494){
10306 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING PY IN Rn07.FOR. A VALUE WILL BE FORCED." << std::endl;
10307 *PY = (AABRA-1.0)*931.4940;
10309 if(std::abs(*PY)>= AABRA*931.494){
10318 std::cout <<
"WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING PZ IN Rn07.FOR. A VALUE WILL BE FORCED." << std::endl;
10319 *PZ = (AABRA-1.0)*931.4940;
10321 if(std::abs(*PZ)>= AABRA*931.494){
10338 v1 = 2.0*
haz(k) - 1.0;
10339 v2 = 2.0*
haz(k) - 1.0;
10340 r = std::pow(v1,2) + std::pow(v2,2);
10343 fac = std::sqrt(-2.*std::log(r)/r);
10345 fgausshaz = v2*fac*sig+xmoy;
10349 fgausshaz=gset*sig+xmoy;
double B(double temperature)
double A(double temperature)
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)
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]