PROGRAM SORT ! ! PROGRAM DESCRIPTION: ! ! This is an example of a program used to sort data ! generated by the statistical model code GEMINI. ! The program determines mass and charge distributions ! of the secondary products. ! ! AUTHORS: ! ! Bob Charity , LBL,GSI,WU ! ! CREATION DATE: 1988 ! ! ! C H A N G E L O G ! ! Date | Name | Description ! ----------------+-------+----------------------------------------------------- ! 19-10-89 | RJC | compound nucleus spin distribution now has a ! diffuseness parameter. see Nucl. Phys. A457 ! (1986)441 ! ----------------+-------+----------------------------------------------------- ! 21-11-89 | RJC | incorporates subroutine Pre_fission into the program ! to calculate pre- and post-fission neutron, proton, ! and alpha-particle multiplicities. This subroutine ! was previously program sort_hinde. ! ----------------+-------+----------------------------------------------------- ! 5-7-90 | RJC | Modified to be compatable with new version of ! GEMINI with new weighting scheme ! ----------------+-------+----------------------------------------------------- ! 22-8-90 | RJC | modified to be compatible with new gemini version ! Also the distribution of pre-fission charge ! particle multiplicities are obtained ! ----------------+-------+----------------------------------------------------- ! [change_entry] ! USE trig USE gem_kind USE sortsub IMPLICIT NONE INTEGER :: L INTEGER :: I_read_error,N_tot,M_frag,i,nn_frag,min_frag INTEGER :: unit_in,unit_out REAL (KIND=r4) :: factor,Jmax,dif REAL (KIND=r4) :: sum_weight,sum_weight_sq REAL (KIND=r4) :: fusion,av,av_error REAL (KIND=r4) :: Z_dist(100), Z_dist_2(100), z_dist_n(100),a_dist(300) REAL (KIND=r4) :: Z_evap(100), Z_evap_2(100), Z_evap_n(100), A_average(100) REAL (KIND=r4) :: A_sig(100), Vz(100), Vzs(100), sumc(100), average(100) REAL (KIND=r4) :: sig(100), spectra_m(15),a_dist_n(300) REAL (KIND=r4) :: Ek(6,0:100),ek2(6,0:100),ang(6,0:180),A_dist_2(300) REAL (KIND=r4) :: sumkp3 , akp3 REAL (KIND=r4) :: kxray REAL (KIND=r4) :: n_residue REAL (KIND=r4) :: charged_pre(0:20),charged_all(0:20) REAL (KIND=r4) :: mean_post_pre(0:20),mean_post_all(0:20) REAL (KIND=r4) :: neut(0:9,3),prot(0:9,3),alph(0:9,3),total(0:9) REAL (KIND=r4) :: deut(0:9,3),trit(0:9,3),he3(0:9,3) REAL (KIND=r4) :: fis_ang(0:36) LOGICAL, SAVE :: first=.true. sumkp3 = 0 akp3=0 N_TOT = 0 M_frag = 0 min_frag = 10000 DO unit_in = 2 unit_out = 8 call start_sort(unit_in,unit_out,I_read_error) if (I_read_error /= 0) EXIT if (first) THEN WRITE (UNIT=*, FMT=*) "pi-lambda-bar" WRITE (UNIT=*, FMT=*) "pi-lambda-bar=0, calculation" READ (UNIT=5, FMT=*) factor IF (factor <= 0.) CALL plb () WRITE (UNIT=6, FMT="(TR2,a16,f7.4)")& "pi-lambda-bar-2=",factor WRITE (UNIT=unit_out, FMT="(TR2,a16,f7.4)")& "pi-lambda-bar-2=",factor factor = factor/para%r_num WRITE (UNIT=*, FMT=*) "Jmax,diffuseness" READ (UNIT=5, FMT=*) JMAX,dif WRITE (UNIT=unit_out, FMT="(TR2,a34,F7.1)")& "cross section calculated for lmax=", Jmax WRITE (UNIT=unit_out, FMT="(TR2,a12,f7.4,TR1,a5)")& "diffuseness=",dif,"h_bar" first = .false. END IF DO call get_event(l,i_read_error) IF (I_read_error /= 0) EXIT IF (L > Jmax + NINT(11.*dif)) THEN EXIT ELSE sum_weight = sum_weight + weight sum_weight_sq = sum_weight_sq + weight**2. N_TOT = N_TOT + 1 NN_frag = N_frag M_frag = MAX(M_frag,NN_frag) min_frag = MIN(min_frag,nn_frag) IF (dif > 0.)& weight = weight *1./(1.+exp((REAL(l)-Jmax)/dif)) fusion = fusion + weight CALL statistics () CALL pre_fission () END IF END DO CLOSE (UNIT=2) END DO WRITE (UNIT=6, FMT="(TR2,a5,i4,a5,TR2,a7,i7,/TR2,a15,i5,A15,I5)")& "lmax=", L, "h-bar", "Events=", N_tot,& "Maximum N_frag=", M_frag,"minimum N_frag=",min_frag WRITE (UNIT=8, FMT="(TR2,a5,i4,a5,TR2,a7,i7,/TR2,a15,i5)")& "lmax=", L, "h-bar", "Events=", N_tot,& "Maximum N_frag=", M_frag CALL output () WRITE (UNIT=6, FMT="(TR2,a24,i7,/TR2,a31,f9.2,TR2,a6,f9.2)")& "Number of Simulations = ", N_tot,& "Sums of Weights of all events =", sum_weight,& "error=",SQRT(sum_weight_sq) WRITE (UNIT=8, FMT="(TR2,a24,i7,/TR2,a31,f9.2,TR2,A6,F9.2)")& "Number of Simulations = ", N_tot,& "Sums of Weights of all events =", sum_weight,& "error=",SQRT(sum_weight_sq) av = sum_weight/REAL(N_tot) av_error = SQRT(sum_weight_sq)/REAL(N_tot) WRITE (UNIT=6, FMT="(TR2,a32,f6.3,TR2,a10,F6.3,/TR2,a19)")& "average weight per event event =", av, "plus-minus",av_error WRITE (UNIT=8, FMT="(TR2,a32,f6.3,TR2,a10,F6.3,/TR2,a19)")& "average weight per event event =", av, "plus-minus",av_error CLOSE (UNIT=8) CONTAINS SUBROUTINE STATISTICS () ! ! FUNCTIONAL DESCRIPTION: ! ! This subroutine increments the histograms for Z-distributions, ! calculates average velocity as function of fragment Z_value ! find coincidences between two heaviest Z"s, average mass as a function ! of Z_value, and multiplicities for many fragment events. ! ! DUMMY ARGUMENTS: ! ! none ! ! IMPLICIT INPUTS: ! ! N_frag - number of fragments in an event ! N_binary - number of binary divisions in a event ! weight - weighting factor for the event ! Z(1:N_frag) - proton number of all fragments in an event ! A(1:N_frag) - nucleon number of all fragments in an event ! v(1:N_frag) - velocity in cm/ns of all fragments in an event ! ! IMPLICIT OUTPUTS: ! ! ak3 - associated multiplicity of complex frag ! sumkp3 - normalization for above ! arrays: ! Z_dist - cross section for each Z species ! Z_dist_2 - error for above ! Z_dist_n - number of simulated events as function of Z-species ! Z_evap - cross section for evaporation residues Z-species ! Z_evap_2 - error for above ! Z_evap_n - simulated events for above ! average - average total charge of coincidence events as function of Z1 ! sig - standard deviation of total charge of coincidence events ! sumc - cross section for coincidence events ! A_average - average A values for each Z-species ! A_sig - stand deviation of A for each Z-species ! Vz - average velocity for each Z_species ! VZs - standard deviation of velocity for each Z-species ! spectra_m - multiplicity spectra ! ! ! SIDE EFFECTS: ! ! [description_or_none] ! ! IMPLICIT NONE INTEGER :: j,izmax1,izmax2,multip_h,mkp3,iz,max1,max2 REAL (KIND=r4) :: ekinetic,zt LOGICAL :: residue iZmax1 = 0. iZmax2 = 0. Multip_h = 0. Mkp3 = 0. residue = .FALSE. IF (frag(n_frag)%id == 1) THEN residue = .TRUE. n_residue = n_residue + weight end If DO I=1,N_FRAG IF (frag(i)%Z == 0 .AND. frag(i)%A > 1)& frag(i)%A = frag(I)%A - 256. !k-xray cross IF (frag(i)%Z == 0 .AND. frag(i)%A == -10) THEN kxray = kxray + weight END IF !energy spectra IF (para%I_angle) THEN IF (frag(i)%Z == 0 .AND. frag(i)%A == 1) THEN ekinetic = 0.5*frag(i)%vel%v**2.*1.042 j = INT(ekinetic) IF (j <= 100) THEN ek(1,j) = ek(1,j) + weight IF (residue) ek2(1,j) = ek2(1,j) + weight END IF j = INT(frag(i)%vel%theta/10) IF (frag(i)%vel%theta <= 1) frag(i)%vel%theta=1 IF (frag(i)%vel%theta >= 179) frag(i)%vel%theta=179 IF (residue) ang(1,j) = ang(1,j)& + weight/SIN_D(frag(i)%vel%theta) ELSE IF (frag(i)%Z == 1 .AND. frag(i)%A == 1) THEN ekinetic = 0.5*frag(i)%vel%v**2.*1.042 j = INT(ekinetic) IF (j <= 100) THEN ek(2,j) = ek(2,j) + weight IF (residue) ek2(2,j) = ek2(2,j) + weight END IF j = INT(frag(i)%vel%theta/10) IF (frag(i)%vel%theta <= 1) frag(i)%vel%theta=1 IF (frag(i)%vel%theta >= 179) frag(i)%vel%theta=179 IF (residue)ang(2,j) = ang(2,j)& + weight/SIN_D(frag(i)%vel%theta) ELSE IF (frag(i)%Z == 2 .AND. frag(i)%A == 4) THEN ekinetic = 0.5*4.*frag(i)%vel%v**2.*1.042 j = INT(ekinetic) IF (j <= 100) THEN ek(3,j) = ek(3,j) + weight IF (residue) ek2(3,j) = ek2(3,j) + weight END IF j = INT(frag(i)%vel%theta/10) IF (frag(i)%vel%theta <= 1) frag(i)%vel%theta=1 IF (frag(i)%vel%theta >= 179) frag(i)%vel%theta=179 IF (residue)ang(3,j) = ang(3,j)& + weight/SIN_D(frag(i)%vel%theta) ELSE IF (frag(i)%Z == 1 .AND. frag(i)%A == 2) THEN ekinetic = 0.5*2.*frag(i)%vel%v**2.*1.042 j = INT(ekinetic) IF (j <= 100) THEN ek(4,j) = ek(4,j) + weight IF (residue) ek2(4,j) = ek2(4,j) + weight END IF j = INT(frag(i)%vel%theta/10) IF (frag(i)%vel%theta <= 1) frag(i)%vel%theta=1 IF (frag(i)%vel%theta >= 179) frag(i)%vel%theta=179 IF (residue) ang(4,j) = ang(4,j)& + weight/SIN_D(frag(i)%vel%theta) ELSE IF (frag(i)%Z == 1 .AND. frag(i)%A == 3) THEN ekinetic = 0.5*3.*frag(i)%vel%v**2.*1.042 j = INT(ekinetic) IF (j <= 100) THEN ek(5,j) = ek(5,j) + weight IF (residue) ek2(5,j) = ek2(5,j) + weight END IF j = INT(frag(i)%vel%theta/10) IF (frag(i)%vel%theta <= 1) frag(i)%vel%theta=1 IF (frag(i)%vel%theta >= 179) frag(i)%vel%theta=179 IF (residue) ang(5,j) = ang(5,j)& + weight/SIN_D(frag(i)%vel%theta) ELSE IF (frag(i)%Z == 2 .AND. frag(i)%A == 3) THEN ekinetic = 0.5*3.*frag(i)%vel%v**2.*1.042 j = INT(ekinetic) IF (j <= 100) THEN ek(6,j) = ek(6,j) + weight IF (residue) ek2(6,j) = ek2(6,j) + weight END IF j = INT(frag(i)%vel%theta/10) IF (frag(i)%vel%theta <= 1) frag(i)%vel%theta=1 IF (frag(i)%vel%theta >= 179) frag(i)%vel%theta=179 IF (residue) ang(6,j) = ang(6,j)& + weight/SIN_D(frag(i)%vel%theta) END IF END IF !calculate multiplicity of heavies IF (frag(i)%Z > 2) Multip_h = Multip_h + 1 if (frag(i)%Z > 2 .AND. frag(i)%Z < 21) Mkp3 = Mkp3 + 1 IF (frag(i)%Z > 0) THEN !increment Charge distributions IZ = frag(I)%Z + 1 Z_DIST(IZ) = Z_DIST(IZ) + weight Z_dist_2(IZ) = Z_dist_2(Iz) + weight**2. Z_dist_n(iz) = Z_dist_n(Iz) + 1 A_average(IZ) = A_average(IZ) + weight*REAL(frag(I)%A) A_sig(IZ) = A_sig(IZ) + weight*REAL(frag(i)%A)**2. Vz(iz) = Vz(iz) + weight*frag(i)%vel%v Vzs(iz) = Vzs(iz) + weight*frag(i)%vel%v**2. a_dist(frag(i)%A) = a_dist(frag(i)%A) + weight a_dist_2(frag(i)%A) = a_dist_2(frag(i)%A) + weight**2. a_dist_n(frag(i)%A) = a_dist_n(frag(i)%A) + 1 ! increment Evaporation charge distributions if (N_binary==0) then Z_EVAP(IZ) = Z_EVAP(IZ) + weight Z_evap_2(IZ) = Z_evap_2(IZ) + weight**2. Z_evap_n(iz) = Z_evap_n(iz) + 1 else ! find two heaviest fragments IF (iZmax1 <= frag(i)%Z) THEN iZmax2 = iZmax1 iZmax1 = frag(i)%Z max2 = max1 max1 = i ELSE IF (iZmax2 < frag(I)%Z) THEN iZmax2 = frag(i)%Z max2 = i END IF END IF END IF END DO !increment coincidence spectra IF (izmax2 >= 2.) THEN sumc(izmax1) = sumc(izmax1) + weight sumc(izmax2) = sumc(izmax2) + weight Zt = REAL(Izmax1+Izmax2) average(izmax1) = average(izmax1) + weight*Zt average(izmax2) = average(izmax2) + weight*Zt sig(izmax1) = sig(izmax1) + weight*Zt**2. sig(izmax2) = sig(izmax2) + weight*Zt**2. END IF !increment multiplicity spectra IF (multip_h > 0) THEN spectra_M(multip_h) = spectra_M(multip_h) + weight END IF IF (mkp3 >= 1) THEN sumkp3 = sumkp3 + weight akp3 = akp3 + weight*REAL(Mkp3) END IF if (frag(n_frag)%z < 50) THEN j = INT(frag(max1)%vel%theta/5.) fis_ang(j) = fis_ang(j) + weight j = INT(frag(max2)%vel%theta/5.) fis_ang(j) = fis_ang(j) + weight END IF RETURN END SUBROUTINE STATISTICS SUBROUTINE OUTPUT () ! ! FUNCTIONAL DESCRIPTION: ! ! This subroutine writes out the contents of the histograms ! ! DUMMY ARGUMENTS: ! ! none ! ! IMPLICIT INPUTS: ! ! factor - normalization factor to give cross sections in mb ! ak3 - associated multiplicity of complex frag ! sumkp3 - normalization for above ! ! arrays: ! Z_dist - cross section for each Z species ! Z_dist_2 - error for above ! Z_dist_n - number of simulated events as function of Z-species ! Z_evap - cross section for evaporation residues Z-species ! Z_evap_2 - error for above ! Z_evap_n - simulated events for above ! average - average total charge of coincidence events as function of Z1 ! sig - standard deviation of total charge of coincidence events ! sumc - cross section for coincidence events ! A_average - average A values for each Z-species ! A_sig - stand deviation of A for each Z-species ! Vz - average velocity for each Z_species ! VZs - standard deviation of velocity for each Z-species ! spectra_m - multiplicity spectra ! neut - pre- and post-fission neutron multiplicity specta ! prot - pre- and post-fission proton multiplicity spectra ! alph - pre- and post-fission alpha-particle multiplicty spectra ! total - normalization for above spectra ! ! IMPLICIT OUTPUTS: ! ! none ! ! ! SIDE EFFECTS: ! ! Writes data to output file ! ! IMPLICIT NONE INTEGER :: k,i,ii,jj,j REAL (KIND=r4) :: mass_asy,dist,dist_2,sumation WRITE (UNIT=6, FMT="(TR2,A17,TR1,F9.2,A2)")& "Fusion Xsection=",fusion*factor,"mb" WRITE (UNIT=8, FMT="(TR2,A17,TR1,F9.2,A2)")& "Fusion Xsection=",fusion*factor,"mb" WRITE (*,*) "fission angular distribution" DO k = 0,35 WRITE (*,*) REAL(k+0.5)*5, fis_ang(k)/factor,fis_ang(k)**0.5/factor END DO k = 100 DO IF (z_dist(k) /= 0) EXIT k = k - 1 END DO WRITE (UNIT=8, FMT="(TR3,a1,TR6,a8,TR7,a8,TR5,a8)")& "Z", "sigma(Z)", "st error", "N events" DO I=1,k IF (Z_dist(i) > 0) THEN dist = z_dist(i)*factor dist_2 = Z_dist_2(i)**0.5*factor WRITE (UNIT=8, FMT="(TR2,I3,TR2,F19.13,TR2,F19.13,TR2,f8.0)")& I-1, DIST, dist_2, Z_dist_n(i) ! WRITE (UNIT=10, FMT="(TR2,I3,A,F14.8)") I-1,",",dist END IF END DO !write out evaporation residue cross sections WRITE (UNIT=8, FMT="(/TR2,a)") "Calssical Evaporation Residues" WRITE (UNIT=8, FMT="(TR3,a1,TR6,a8,TR7,a8,TR5,a8)")& "Z", "sigma(Z)", "st error", "N events" DO i=5,k IF (Z_evap(i) > 0) THEN Z_evap(i) = Z_evap(i)*factor Z_evap_2(i) = Z_evap_2(I)**0.5*factor WRITE (UNIT=8, FMT="(TR2,I3,TR2,F12.6,TR2,F12.6,TR2,f12.4)")& I-1, Z_evap(I), Z_evap_2(i), Z_evap_n(i) ! WRITE (UNIT=10, FMT="(TR2,I3,A,F14.8)") I-1,",",Z_evap(i) END IF END DO k = 300 DO IF (A_dist(k) == 0) EXIT k = k - 1 END DO WRITE (UNIT=8, FMT="(TR3,a1,TR6,a8,TR7,a8,TR5,a8)")& "A", "sigma(A)", "st error", "N events" DO I=1,k IF (A_dist(i) > 0) THEN dist = A_dist(i)*factor dist_2 = A_dist_2(i)**0.5*factor WRITE (UNIT=8, FMT="(TR2,I3,TR2,F19.13,TR2,F19.13,TR2,f8.0)")& I, DIST, dist_2, A_dist_n(i) ! WRITE (UNIT=10, FMT="(TR2,I3,A,F14.8)") I-1,",",dist END IF END DO !write coincidence results WRITE (UNIT=8, FMT="(/TR2,a,/TR2,a,/TR2,a,TR3,a1,TR6,a7,TR7,a12)")& "For coincidences between the 2 heaviest fragments.",& "Charge lost due to evaporation before and after scission",& "as a function of the Z of one fragment",& "Z", "", "sigma(Zloss)" ii=5 DO IF (sumc(ii) <= 0.) EXIT average(ii) = average(ii)/sumc(ii) sig(ii) = sig(ii)/sumc(ii) sig(ii) = sig(ii) - average(ii)**2. if (sig(ii) <= 0.) sig(ii) = 0. sig(ii) = sig(ii)**0.5 average(ii) = para%Zcn - Average(ii) WRITE (UNIT=8, FMT="(TR2,I3,TR7,f7.3,TR10,f5.3)")& ii, average(ii), sig(ii) ii = ii + 1 IF (sumc(ii) == 0) ii = ii + 1 END DO ! WRITE (UNIT=8, FMT="(/TR2,a,TR3,a1,TR7,a3,TR6,a6,TR4,a21)")& ! "Average mass of each Z-species",& ! "Z", "", "sig(A)", "-2.08Z+0.0029*Z**2" ! DO I=1,k ! IF (Z_dist(i) > 0) THEN ! a_average(i) = a_average(i)/Z_dist(i) ! a_sig(i) = a_sig(i)/Z_dist(i) ! a_sig(i) = a_sig(i) - a_average(i)**2. ! if (A_sig(i)%le.0.) a_sig(i) = 0. ! a_sig(i) = a_sig(i)**0.5 ! da = a_average(i) - 2.08*REAL(i-1)& ! -0.0029*REAL(i-1)**2. ! WRITE (UNIT=8, FMT="(TR2,I3,TR4,F6.2,TR4,F6.2,TR11,f6.2)")& ! i-1, a_average(i), a_sig(i), da ! END IF ! END DO IF (para%I_angle) THEN WRITE (UNIT=8, FMT="(/TR2,a,/TR3,a1,TR7,a5,TR6,a8)")& "Average c-m velocity of each Z-species",& "Z", "", "sig(Vcm)" DO I=4,k IF (Z_dist(i) > 0) THEN Vz(i) = Vz(i)/Z_dist(i) Vzs(i) = Vzs(i)/Z_dist(i) Vzs(i) = Vzs(i) - Vz(i)**2. if (Vzs(i) <= 0.) Vzs(i) = 0. Vzs(i) = Vzs(i)**0.5 WRITE (UNIT=8, FMT="(TR2,I3,TR4,F6.2,TR4,F6.2)")& i-1, Vz(i), Vzs(i) END IF END DO END IF WRITE (UNIT=8, FMT="(/TR2,a,/TR2,a12,TR8,a11)")& "Cross sections for heavy particle Multiplicities",& "Multiplicity", "xsection mb" DO i=1,10 WRITE (UNIT=8, FMT="(TR7,I2,TR16,E10.4)")& I, spectra_m(i)*factor END DO IF (sumkp3 > 0.) THEN akp3 = akp3/sumkp3 WRITE (UNIT=8, FMT="(/TR2,a,/TR2,a,f5.2)")& "average multiplicity of IMFs (2 0.) THEN DO jj=1,3 neut(i,jj) = neut(i,jj)/total(i) prot(i,jj) = prot(i,jj)/total(i) alph(i,jj) = alph(i,jj)/total(i) deut(i,jj) = deut(i,jj)/total(i) trit(i,jj) = trit(i,jj)/total(i) he3(i,jj) = he3(i,jj)/total(i) END DO mass_asy = REAL(i)/10. WRITE (UNIT=8, FMT="(TR2, f4.2, 3(TR2, 3(TR1, f6.3)))")& mass_asy, neut(i,1), prot(i,1), alph(i,1),& neut(i,2), prot(i,2), alph(i,2),& neut(i,3), prot(i,3), alph(i,3) WRITE (UNIT=8, FMT="(TR2,a3,3(TR2,3(TR4,a3)))")& "asy","d","t","He3","d","t","He3","d","t","He3" WRITE (UNIT=8, FMT="(TR2, f4.2, 3(TR2, 3(TR1, f6.3)))")& mass_asy, deut(i,1), trit(i,1), he3(i,1),& deut(i,2), trit(i,2), he3(i,2),& deut(i,3), trit(i,3), he3(i,3) END IF END DO sumation = 0. DO i=0,10 sumation = sumation + charged_pre(i) END DO WRITE (UNIT=8, FMT="(TR2,a,TR2,a,TR2,a,TR2,a)")& "charged particles","pre-fission","all","mean_post" DO i=0,20 IF (charged_pre(i) > 0.) THEN mean_post_pre(i) = mean_post_pre(i)/charged_pre(i) ELSE mean_post_pre(i) = 0. END IF IF (charged_all(i) > 0.) THEN mean_post_all(i) = mean_post_all(i)/charged_all(i) ELSE mean_post_all(i) = 0. END IF IF (sumation > 0) WRITE (UNIT=8,& FMT="(TR2,I2,TR2,F8.3,TR2,F8.3,TR2,F8.3,TR2,F8.3)")& i,charged_pre(i)/sumation,charged_all(i)/sumation,& mean_post_pre(i),mean_post_all(i) END DO WRITE (UNIT=8, FMT="(TR2,A24,TR2,E15.5,TR1,A2)")& "k xray-fission xsection=",kxray*factor,"mb" IF (para%i_angle) THEN WRITE (UNIT=8, FMT=*) WRITE (unit=8, FMT="(TR2,A)")& "kinetic energy, d_sigma/de for n,p,a" DO j=0,100 WRITE (UNIT=8, FMT="(TR2,F5.1,3(TR2,E10.5))")& REAL(j)+.5,ek(1,j)*factor,ek(2,j)*factor,ek(3,j)*factor END DO WRITE (UNIT=8, FMT=*) WRITE (unit=8, FMT="(TR2,A)")& "kinetic energy, d_sigma/de for d,t,3he" DO j=0,100 WRITE (UNIT=8, FMT="(TR2,F5.1,3(TR2,E10.5))")& REAL(j)+.5,ek(4,j)*factor,ek(5,j)*factor,ek(6,j)*factor END DO WRITE (UNIT=8, FMT=*) WRITE (unit=8, FMT="(TR2,A)")& "kinetic energy, d_sigma/de for n,p,a from residue" DO j=0,100 WRITE (UNIT=8, FMT="(TR2,F5.1,3(TR2,E10.5))")& REAL(j)+.5,ek2(1,j)/n_residue/12.56,& ek2(2,j)/n_residue/12.56,& ek2(3,j)/n_residue/12.56 END DO WRITE (UNIT=8, FMT=*) WRITE (unit=8, FMT="(TR2,A)")& "kinetic energy, d_sigma/de for d,t,3he from residue" DO j=0,100 WRITE (UNIT=8, FMT="(TR2,F5.1,3(TR2,E10.5))")& REAL(j)+.5,ek2(4,j)/n_residue/12.56,& ek2(5,j)/n_residue/12.56,& ek2(6,j)/n_residue/12.56 END DO WRITE (UNIT=8, FMT=*) WRITE (unit=8, FMT="(TR2,A)")& "angular distribtion relatize to spin axis n,p,a" DO j=0,17 WRITE (UNIT=8, FMT="(TR2,I3,3(TR2,E10.5))")& j*10+5,ang(1,j)/10*180/3.14159/6.2832/n_residue& ,ang(2,j)/10*180/3.14159/6.2832/n_residue& ,ang(3,j)/10*180/3.14159/6.2832/n_residue END DO WRITE (UNIT=8, FMT=*) WRITE (unit=8, FMT="(TR2,A)")& "angular distribtion relatize to spin axis d,t,3he" DO j=0,17 WRITE (UNIT=8, FMT="(TR2,I3,3(TR2,E10.5))")& j*10+5,ang(4,j)/10*180/3.14159/6.2832/n_residue& ,ang(5,j)/10*180/3.14159/6.2832/n_residue& ,ang(6,j)/10*180/3.14159/6.2832/n_residue END DO END IF RETURN END SUBROUTINE OUTPUT SUBROUTINE plb () ! ! FUNCTIONAL DESCRIPTION: ! ! calculates the quantity pi-lambda-bar squared for use in normalising ! calculated yields into cross sections ! ! DUMMY ARGUMENTS: ! ! output: factor (REAL (KIND=r4) ::) pi-lambda-bar-2 in mb ! ! IMPLICIT INPUTS: ! ! none ! ! IMPLICIT OUTPUTS: ! ! none ! ! ! SIDE EFFECTS: ! ! none ! ! IMPLICIT NONE REAL (KIND=r4) :: Ap,AT,EPA,U !read in projectile and target mass and beam energy WRITE (UNIT=*, FMT=*) "Ap,At,E/A" READ (UNIT=5, FMT=*) Ap,At,EPA u = Ap*At/(Ap+At) factor = 656.80/(U**2.*Epa) RETURN END SUBROUTINE plb SUBROUTINE pre_fission () ! ! FUNCTIONAL DESCRIPTION: ! ! This subroutine gives the ! pre- and post fission neutron multiplicities for neutrons, protrons, ! and alpha particles. For each event, the program first finds the ! heaviest two fragments. (purely evaporative events are ignored)% No ! check is made to see if these are "fission fragments", i.e. they may ! not be emitted 180 degress apart in the centre-of-mass. The later may ! be of importance for the higher excitation energies where ternary and ! higher order events have significant cross sections. However, for the ! lower excitation energies, these fragment will be fission-fragments. ! Information on the decay history of these fragments is contained in the ! integer variable I_decay. However it is not so readily understandable. ! This program constructs the array hist(i,j) [i=fragment number and j= ! family tree level) for each fragment giving its ancestors, in order, in ! the family tree or decay chain. By comparing the Hist array for the two ! fission fragments we can see their common ancestor, the fissioning ! system. All fragments (light or otherwise) emitted from this fissioning ! system or its ancestors are considered pre-fission. And simularly the ! other fragments can be classified as post-fission from either the ! light or heavier fragment. The total mass of all the post-fission ! particles is added up to give the primary fission-fragments. - this ! may lead to a redefinition of which is the heavy and light ! fission-fragment. The primary mass asymmetry of the fission fragments ! (A1-A2)/(A1+A2) is calculated and the multiplicities are sorted ! acccordingly. The multiplicities for heavier fragemnts has not been ! sorted but could easily be incorporated in this program ! ! DUMMY ARGUMENTS: ! ! none ! ! IMPLICIT INPUTS: ! ! N_frag - number of fragments in an event ! N_binary - number of binary divisions in a event ! weight - weighting factor for the event ! Z(1:N_frag) - proton number of all fragments in an event ! A(1:N_frag) - nucleon number of all fragments in an event ! v(1:N_frag) - velocity in cm/ns of all fragments in an event ! I_decay(1:N_frag) - branch of decay tree for each fragment ! ! IMPLICIT OUTPUTS: ! ! neut - pre- and post-fission neutron multiplicity specta ! prot - pre- and post-fission proton multiplicity spectra ! alph - pre- and post-fission alpha-particle multiplicty spectra ! total - normalization for above spectra ! ! ! SIDE EFFECTS: ! ! none ! ! ! C H A N G E L O G ! ! Date | Name | Description ! ----------------+-------+----------------------------------------------------- ! 16-11-89 | RJC | the logical conditions for determining the the HIST ! arrays were not strictly correct. Some many-fragment ! events were not treated correctly. Also some ! beutification was made (lables removed) ! ----------------+-------+----------------------------------------------------- IMPLICIT NONE INTEGER :: hist(50,8) INTEGER :: A_add(3),iamax1,iamax2,i,k,jj,kk,i1,i2 INTEGER :: n_mul(3),p_mul(3),a_mul(3),d_mul(3),t_mul(3),he3_mul(3) REAL (KIND=r4) :: mass_asy INTEGER :: pre,all iamax1 = 0. iamax2 = 0. IF (N_binary > 0) THEN DO i=1,N_frag IF (frag(i)%Z < 0) CYCLE DO k = 1,8 hist(i,k) = 0 END DO jj = 8 DO IF (IBITS (frag(i)%id, jj-1, 1) /= 0) EXIT jj = jj - 1 END DO kk = 0 ! determine history of fragment, i.e. from what intermediate fragments was it ! emitted DO k=jj,1 ,-1 kk = kk + 1 hist(i,kk) = IBITS (frag(i)%id, k-1, 9-k) END DO IF (iamax1 <= frag(i)%A) THEN iamax2 = iamax1 iamax1 = frag(i)%A i2 = i1 i1 = i ELSE IF (iamax2 < frag(i)%A) THEN iamax2 = frag(i)%A i2 = i END IF END DO IF (iamax2 > 4) THEN ! find where histories diverge I = 1 DO IF (hist(I1,I) /= hist(I2,I)) EXIT I = I + 1 END DO A_add = 0 n_mul = 0 p_mul = 0 a_mul = 0 d_mul = 0 t_mul = 0 he3_mul = 0 DO jj=1,N_frag if (frag(jj)%Z < 0) CYCLE !ignore gammas and xrays IF (hist(jj,i) == hist(i1,i)) THEN k = 2 ! post fission from heavy fragment ELSE IF (hist(jj,i) == hist(i2,i)) THEN k = 3 ! post fission from light fragment ELSE k = 1 !pre-fission END IF IF (k /= 4) THEN SELECT CASE(frag(jj)%z) CASE (0) !neutron if (frag(jj)%A == 1) n_mul(k) = n_mul(k) + 1 CASE (1) ! H SELECT CASE (frag(jj)%A) CASE (1) !proton p_mul(k) = p_mul(k) + 1 CASE (2) !deuteron d_mul(k) = d_mul(k) + 1 CASE (3) !triton t_mul(k) = t_mul(k) + 1 END SELECT CASE (2) !He SELECT CASE (frag(jj)%A) CASE (3) he3_mul(k) = he3_mul(k) + 1 CASE (4) !alpha particle a_mul(k) = a_mul(k) + 1 END SELECT END SELECT END IF A_add(k) = A_add(k) + MAX(0,frag(jj)%A) END DO IF (A_add(2)+A_add(3) > 0) THEN mass_asy = ABS(REAL(A_add(2)-A_add(3))& /REAL(A_add(2)+A_add(3))) k = NINT(mass_asy*10.) neut(k,1) = neut(k,1) + REAL(n_mul(1))*weight prot(k,1) = prot(k,1) + REAL(p_mul(1))*weight alph(k,1) = alph(k,1) + REAL(a_mul(1))*weight deut(k,1) = deut(k,1) + REAL(d_mul(1))*weight trit(k,1) = trit(k,1) + REAL(t_mul(1))*weight he3(k,1) = he3(k,1) + REAL(he3_mul(1))*weight total(k) = total(k) + weight IF (A_add(3) <= A_add(2)) THEN i1 = 2 i2 = 3 ELSE i1 = 3 i2 = 2 END IF neut(k,2) = neut(k,2) + REAL(n_mul(i1))*weight prot(k,2) = prot(k,2) + REAL(p_mul(i1))*weight alph(k,2) = alph(k,2) + REAL(a_mul(i1))*weight deut(k,2) = deut(k,2) + REAL(d_mul(i1))*weight trit(k,2) = trit(k,2) + REAL(t_mul(i1))*weight he3(k,2) = he3(k,2) + REAL(he3_mul(i1))*weight neut(k,3) = neut(k,3) + REAL(n_mul(i2))*weight prot(k,3) = prot(k,3) + REAL(p_mul(i2))*weight alph(k,3) = alph(k,3) + REAL(a_mul(i2))*weight deut(k,3) = deut(k,3) + REAL(d_mul(i2))*weight trit(k,3) = trit(k,3) + REAL(t_mul(i2))*weight he3(k,3) = he3(k,3) + REAL(he3_mul(i2))*weight IF (mass_asy < 0.71) THEN pre = p_mul(1) + a_mul(1) charged_Pre(pre) = charged_pre(pre) + weight all = pre + p_mul(2) + a_mul(2) + p_mul(3) + a_mul(3) charged_all(all) = charged_all(all) + weight mean_post_pre(pre) = mean_post_pre(pre)& + weight*REAL(all-pre) mean_post_all(all) = mean_post_all(all)& + weight*REAL(all-pre) END IF END IF END IF END IF RETURN END SUBROUTINE pre_fission END PROGRAM SORT