PROGRAM sort_er
USe gem_kind
USE sortsub
USE functs
USE trig
IMPLICIT NONE


REAL (KIND=r4) :: Z_dist(100), Z_dist_2(100), z_dist_n(100)
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)
REAL (KIND=r4) :: Ek(14,0:100),ek2(14,0:100),ang(14,0:180),ang_res(0:50,-1:5)
REAL (KIND=r4) :: ek3(14,0:100),ek4(14,0:100),xsec_er(-1:5)
REAL (KIND=r4) :: ang_res_step(7),A_dist(200),phi_aa(0:90)
REAL (KIND=r4) :: fusion,sumzer,sumner
REAL (KIND=r4) :: kxray,ran,sum_gam_energy=0,sum_missing=0

!prefission mult
!REAL (KIND=r4), DIMENSION(0:20) :: charged_pre,charged_all
!REAL (KIND=r4), DIMENSION(0:20) :: mean_post_pre,mean_post_all
!REAL (KIND=r4), DIMENSION(0:9,3) :: neut,prot,alph
!REAL (KIND=r4), DIMENSION(0:9)   :: total


INTEGER :: L
INTEGER :: I_read_error,N_tot,M_frag,nn_frag,min_frag
INTEGER :: unit_in,unit_out,seed
REAL (KIND=r4) :: factor,Jmax,dif,jmin
REAL (KIND=r4) :: sum_weight,sum_weight_sq
REAL (KIND=r4) :: av,av_error,n_ins,n_all

REAL (KIND=r4) ::  akp3=0.
REAL (KIND=r4) ::  sumkp3=0.
REAL (KIND=r4) :: vs,theta_lab_max,theta_lab_min,z_ev,theta_scatter
REAL (KIND=r4) :: n_residue,mult_n,mult_p,mult_a,mult_d,mult_t
REAL (KIND=r4) :: ek_n,ek_p,ek_a,ek_d,ek_t,mult_3he,ek_3he,mult_6li,mult_7li
REAL (KIND=r4) :: ek_6li,ek_7li,ek2_p,ek2_n,ek2_a,ek2_d,ek2_t,ek2_3he
REAL (KIND=r4) :: ek2_6li,ek2_7li,mult_7be,mult_6he,ek_7be,ek2_7be,ek_6he
REAL (KIND=r4) :: mult_9be,ek_9be,ek2_9be,mult_8li,ek_8li,ek2_8li,ek2_6he
REAL (KIND=r4) :: mult_10be,ek_10be,ek2_10be,mult_8be,ek_8be,ek2_8be
REAL (KIND=r4), DIMENSION(0:120) :: v_residue,av_residue
REAL (KIND=r4), DIMENSION(220) :: a_residue
LOGICAL, SAVE :: first=.true.


seed = 1312635



N_TOT = 0
M_frag = 0
min_frag = 10000
sum_weight = 0.
sum_weight_sq = 0.
fusion = 0.
n_all = 0.
ang_res_step = 0.
n_residue = 0.
mult_n = 0.
mult_p = 0.
mult_d = 0.
mult_t = 0.
mult_3he = 0.
mult_a = 0.
mult_6he = 0.
mult_6li = 0.
mult_7li = 0.
mult_8li = 0.
mult_8be = 0.
mult_9be = 0.
mult_10be = 0.
ek_n = 0.
ek_p = 0.
ek_d = 0.
ek_t = 0.
ek_3he = 0.
ek_a = 0.
ek_6he = 0.
ek_6li = 0.
ek_7li = 0.
ek_8li = 0.
ek_8Be = 0.
ek_9be = 0.
ek_10be = 0.
a_residue = 0.
v_residue = 0.
av_residue = 0.
z_dist = 0.
z_dist_2 = 0.
z_dist_n = 0.
z_evap = 0.
z_evap_2 = 0.
z_evap_n = 0
A_average = 0.
a_sig = 0.
vz = 0.
vzs = 0.
sumc = 0.
average = 0.
sig = 0.
spectra_m = 0.
ek = 0.
ek2 = 0.
ang = 0.
ang_res = 0.
ek3 = 0.
ek4 = 0.
xsec_er = 0.
A_dist = 0.
phi_aa = 0.
sumzer = 0.
sumner = 0.
kxray = 0.
ek2_p = 0.
ek2_n = 0.
ek2_a = 0.
ek2_d = 0.
ek2_t = 0.
ek2_3he = 0.
ek2_6he = 0.
ek2_6li = 0.
ek2_7li = 0.
ek2_8li = 0.
ek2_7be = 0.
ek2_8be = 0.
ek2_9be = 0.
ek2_10be = 0.
n_ins = 0

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=*)&
        "vs,theta_lab_max,theta_lab_min,z_ev,theta_scatter"
        READ (UNIT=5, FMT=*) vs,theta_lab_max,theta_lab_min,z_ev,theta_scatter

        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

        WRITE (UNIT=*, FMT=*) "Jmax,diffuseness,jmin"
        READ (UNIT=5, FMT=*) JMAX,dif,jmin
        WRITE (UNIT=unit_out, FMT="(TR2,a34,F7.1)")&
        "cross section calculated for lmax=", Jmax
        WRITE (UNIT=unit_out, FMT="(TR2,a12,f7.4,1x,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
         if (l < jmin) weight=0
         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.
!
! 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

INTEGER :: j,alpha,i,k,izmax1,izmax2,multip_h,mkp3,ii,iz,ix,jj
REAL (KIND=r4) :: THETA_LAB_Z,zt
REAL (KIND=r4) :: theta_Lab,vxx,vzz,vv
REAL (KIND=r4) :: vx1,vy1,vz1,vx2,vy2,vz2,delta_v,e_rel,vyy,theta

REAL (KIND=r4) :: ek_h(2,100,9)
INTEGER :: mult_aa
REAL (KIND=r4) :: phi_a(30),delta_phi

REAL (KIND=r4) :: d_theta,d_phi,theta_after,phi_after,ekinetic

logical :: residue,alone,primary

iZmax1 = 0
iZmax2 = 0
Multip_h = 0
Mkp3 = 0





residue = .FALSE.
IF (frag(n_frag)%Z >= z_ev) THEN
   n_all = n_all + weight
   vzz = frag(n_frag)%vel%v*COS_D(frag(n_frag)%vel%theta) + vs
   vxx = frag(n_frag)%vel%v*SIN_D(frag(n_frag)%vel%theta)
   vv = (vxx**2.+vzz**2.)**0.5
   if (vv > 0) then
       theta_lab = ACOS_D(vzz/vv)
   else
       theta_lab = 0
       WRITE (UNIT=*, FMT=*) "res"
   END IF


   IF (theta_scatter > 0) THEN
       CALL RANDOM_NUMBER(ran)
       d_theta = 1.128*theta_scatter*(-LOG(ran))**0.5
       CALL RANDOM_NUMBER(ran)
       d_phi = 360.*RAN
       CALL transform(theta_lab,frag(n_frag)%vel%phi,&
       d_theta,d_phi,theta_after,phi_after)

       theta_lab = theta_after
   END IF


   j = theta_lab
   IF (j <= 50) THEN
       IF (theta_lab <= .05) theta_lab = .05
!              ang_res(j,1) = ang_res(j,1) + weight/sin_d(theta_lab)
!               xsec_er(1) = xsec_er(1) + weight
       IF (theta_lab <= 1.47) THEN
           ang_res_step(1) = ang_res_step(1) + weight
       ELSE IF (theta_lab <= 2.15) THEN
           ang_res_step(2) = ang_res_step(2) + weight
       ELSE IF (theta_lab <= 2.97) THEN
           ang_res_step(3) = ang_res_step(3) + weight
       ELSE IF (theta_lab <= 3.84) THEN
           ang_res_step(4) = ang_res_step(4) + weight
       ELSE IF (theta_lab <= 4.81) THEN
           ang_res_step(5) = ang_res_step(5) + weight
       ELSE IF (theta_lab <= 5.76) THEN
           ang_res_step(6) = ang_res_step(6) + weight
       ELSE IF (theta_Lab <= 6.73) THEN
           ang_res_step(7) = ang_res_step(7) + weight
       END IF
   END IF


   if (theta_lab  <= theta_lab_min) n_ins = n_ins + weight


   IF (theta_lab <= theta_lab_max&
   .and. theta_lab >= theta_lab_min) THEN
       residue = .TRUE.
       N_residue = n_residue + weight
       sum_gam_energy = sum_gam_energy + weight*frag(n_frag)%exx
       sum_missing = sum_missing + weight*frag(n_frag)%exx + &
       0.5*REAL(frag(n_frag)%a)*frag(n_frag)%vel%v**2*1.042
       sumzer = sumzer + REAL(frag(n_frag)%z)*weight
       sumner = sumner + REAL(frag(n_frag)%a-frag(n_frag)%z)*weight
       k = INT(vv*40)
       if (k <= 120)THEN
           V_residue(k) = v_residue(k) + weight
           av_residue(k) = av_residue(k) + weight*REAL(frag(n_frag)%a)
       END IF
       A_residue(frag(n_frag)%a) = A_residue(frag(n_frag)%a) + weight
   END IF
END IF

 mult_aa = 0
 ALPHA = 0

DO I=1,n_frag

!k-xray cross
 IF (frag(i)%Z == -1 .AND. frag(i)%a > 10) THEN
     kxray = kxray + weight
 ELSE IF (frag(i)%Z == -1 .AND. frag(i)%A <= 2)THEN !gamma ray
     sum_gam_energy = sum_gam_energy + weight*frag(i)%exx
     sum_missing = sum_missing + weight*frag(i)%exx
   !
 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
         if (residue) THEN
              mult_n = mult_n + weight
              ek_n = ek_n + ekinetic*weight
              ek2_n = ek2_n + ekinetic**2.*weight
         end if
         j = INT(ekinetic)
         IF (j <= 100) THEN
             ek(1,j) = ek(1,j) + weight
             IF (residue) ek2(1,j) = ek2(1,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 0)&
             ek3(1,j) = ek3(1,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 1)&
             ek4(1,j) = ek4(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
         if (residue)then
             mult_p = mult_p + weight
             EK_P = EK_P + ekinetic*weight
             ek2_p = ek2_p + ekinetic**2.*weight
         END IF
         j = INT(ekinetic)
         IF (j <= 100) THEN
             ek(2,j) = ek(2,j) + weight
             IF (residue) ek2(2,j) = ek2(2,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 0)&
             ek3(2,j) = ek3(2,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 1)&
             ek4(2,j) = ek4(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
         ALPHA = alpha + 1
         ekinetic = 0.5*4.*frag(i)%vel%v**2.*1.042

!see if alpha is from a he5, if so decrease energy by 1 mev
!to account for coulomb effect
         if (i > 1 .AND. MODULO(frag(i)%id,2) == 0) THEN
             if (frag(i-1)%a == 1 .AND. frag(i-1)%z == 0&
             .and. frag(i-1)%id == frag(i)%id) THEN
               alone = .true.
               DO ii=1,i-2
                 if (frag(ii)%id  == frag(i)%id) THEN
                    alone = .FALSE.
                    EXIT
                 END IF
               END DO
               if (alone) ekinetic = ekinetic + 1.
             END IF
         END IF
         if (residue) THEN
             mult_a = mult_a + weight
             ek_a = ek_a + ekinetic*weight
             ek2_a = ek2_a + ekinetic**2.*weight
         END IF
         if (residue .AND. frag(i)%vel%theta > 45 .and.&
         frag(i)%vel%theta < 135)then
             mult_aa = mult_aa + 1
             phi_a(mult_aa) = frag(i)%vel%phi
         END IF
         j = INT(ekinetic)
         IF (j <= 100) THEN
             ek(3,j) = ek(3,j) + weight

             IF (residue ) ek2(3,j) = ek2(3,j) + weight

             if (residue .AND. MODULO(frag(i)%id,2) == 0)&
             ek3(3,j) = ek3(3,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 1)&
             ek4(3,j) = ek4(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)

         if (i > 1 .AND. MODULO(frag(i)%id,2) == 0) THEN
             if (frag(i-1)%Z == 2 .AND. frag(i-1)%A == 4 .AND. &
             frag(i-1)%id == frag(i)%id) THEN

             vx1 = frag(i)%vel%v*SIN_D(frag(i)%vel%theta) &
             *COS_D(frag(i)%vel%phi)
             vx2 = frag(i-1)%vel%v*SIN_D(frag(i-1)%vel%theta)&
             *COS_D(frag(i-1)%vel%phi)
             vy1 = frag(i)%vel%v*SIN_D(frag(i)%vel%theta) &
             *SIN_D(frag(i)%vel%phi)
             vy2 = frag(i-1)%vel%v*SIN_D(frag(i-1)%vel%theta)&
             *SIN_D(frag(i-1)%vel%phi)
             vz1 = frag(i)%vel%v*COS_D(frag(i)%vel%theta)
             vz2 = frag(i-1)%vel%v*COS_D(frag(i-1)%vel%theta)

             delta_v = (vx1-vx2)**2. + (vy1-vy2)**2. + (vz1-vz2)**2.

             e_rel = delta_v/.9784**2.*2./2.

             if (e_rel <= .15 .AND. e_rel > 0) THEN
                 primary = .TRUE.

                 DO jj=1,i-2
                     if (frag(jj)%id == frag(i)%id) THEN
                         primary = .FALSE.
                         EXIT
                     END IF
                 END DO

                 vxx = (vx1+vx2)/2.
                 vyy = (vy1+vy2)/2.
                 vzz = (vz1+vz2)/2.
                 vv = (vxx**2.+vyy**2.+vzz**2.)**0.5
                 ekinetic = 0.5*8.*vv**2.*1.042
                 theta = ACOS_D(vzz/vv)
                 if (residue)THEN
                      mult_8be = mult_8be + weight
                      ek_8be = ek_8be + ekinetic*weight
                      ek2_8be = ek2_8be + ekinetic**2.*weight
                 END IF
                 j = INT(ekinetic)
                 IF (j <= 100) THEN
                 ek(14,j) = ek(14,j) + weight
                 IF (residue) ek2(14,j) = ek2(14,j) + weight
                 if (residue .AND. .not. primary)&
                 ek3(14,j) = ek3(14,j) + weight
                 if (residue .AND. primary)&
                 ek4(14,j) = ek4(14,j) + weight
              END IF
              j = INT(theta/10)
              IF (theta <= 1) theta=1
              IF (theta >= 179)theta=179
              IF (residue)ang(14,j) = ang(14,j)&
              + weight/sin_d(theta)

              END IF

            END IF
         END IF
     ELSE IF (frag(i)%Z == 1 .AND. frag(i)%a == 2) THEN
         ekinetic = 0.5*2.*frag(i)%vel%v**2.*1.042
         if (residue)THEN
              mult_d = mult_d + weight
              ek_d = ek_d + ekinetic*weight
              ek2_d = ek2_d + ekinetic**2.*weight
         END IF
         j = INT(ekinetic)
         IF (j <= 100) THEN
             ek(4,j) = ek(4,j) + weight
             IF (residue) ek2(4,j) = ek2(4,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 0)&
             ek3(4,j) = ek3(4,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 1)&
             ek4(4,j) = ek4(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
         if (residue)THEN
             mult_t = mult_t + weight
             ek_t = ek_t + ekinetic*weight
             ek2_t = ek2_t + ekinetic**2.*weight
         END IF
         j = INT(ekinetic)
         IF (j <= 100) THEN
             ek(5,j) = ek(5,j) + weight
             IF (residue) ek2(5,j) = ek2(5,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 0)&
             ek3(5,j) = ek3(5,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 1)&
             ek4(5,j) = ek4(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 (residue)THEN
             mult_3he = mult_3he + weight
             ek_3he = ek_3he + ekinetic*weight
             ek2_3he = ek2_3he + ekinetic**2.*weight
         END IF

         IF (j <= 100) THEN
             ek(6,j) = ek(6,j) + weight
             IF (residue) ek2(6,j) = ek2(6,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 0)&
             ek3(6,j) = ek3(6,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 1)&
             ek4(6,j) = ek4(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)
     ELSE IF (frag(i)%Z == 3 .AND. frag(i)%a == 6) THEN
         ekinetic = 0.5*6.*frag(i)%vel%v**2.*1.042
         j = INT(ekinetic)
         if (residue)THEN
             mult_6li = mult_6li + weight
             ek_6li = ek_6li + ekinetic*weight
             ek2_6li = ek2_6li + ekinetic**2.*weight
         END IF
         IF (j <= 100) THEN
             ek(7,j) = ek(7,j) + weight
             IF (residue) ek2(7,j) = ek2(7,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 0)&
             ek3(7,j) = ek3(7,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 1)&
             ek4(7,j) = ek4(7,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(7,j) = ang(7,j)&
         + weight/SIN_D(frag(i)%vel%theta)

     ELSE IF (frag(i)%Z == 3 .AND. frag(i)%a == 7) THEN
         ekinetic = 0.5*7.*frag(i)%vel%v**2.*1.042
         j = INT(ekinetic)
         if (residue)THEN
             mult_7li = mult_7li + weight
             ek_7li = ek_7li + ekinetic*weight
             ek2_7li = ek2_7li + ekinetic**2.*weight
         END IF
         IF (j <= 100) THEN
             ek(8,j) = ek(8,j) + weight
             IF (residue) ek2(8,j) = ek2(8,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 0)&
             ek3(8,j) = ek3(8,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 1)&
             ek4(8,j) = ek4(8,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(8,j) = ang(8,j)&
         + weight/SIN_D(frag(i)%vel%theta)

     ELSE IF (frag(i)%Z == 3 .AND. frag(i)%a == 8) THEN
         ekinetic = 0.5*8.*frag(i)%vel%v**2.*1.042
         j = INT(ekinetic)
         if (residue)THEN
             mult_8li = mult_8li + weight
             ek_8li = ek_8li + ekinetic*weight
             ek2_8li = ek2_8li + ekinetic**2.*weight
         END IF
         IF (j <= 100) THEN
             ek(9,j) = ek(9,j) + weight
             IF (residue) ek2(9,j) = ek2(9,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 0)&
             ek3(9,j) = ek3(9,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 1)&
             ek4(9,j) = ek4(9,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(9,j) = ang(9,j)&
         + weight/SIN_D(frag(i)%vel%theta)

     ELSE IF (frag(i)%Z == 4 .AND. frag(i)%a == 7) THEN
         ekinetic = 0.5*7.*frag(i)%vel%v**2.*1.042
         j = INT(ekinetic)
         if (residue)THEN
             mult_7Be = mult_7be + weight
             ek_7be = ek_7be + ekinetic*weight
             ek2_7be = ek2_7be + ekinetic**2.*weight
         END IF
         IF (j <= 100) THEN
             ek(10,j) = ek(10,j) + weight
             IF (residue) ek2(10,j) = ek2(10,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 0)&
             ek3(10,j) = ek3(10,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 1)&
             ek4(10,j) = ek4(10,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(10,j) = ang(10,j)&
         + weight/SIN_D(frag(i)%vel%theta)

     ELSE IF (frag(i)%Z == 4 .AND. frag(i)%a == 10) THEN
         ekinetic = 0.5*10.*frag(i)%vel%v**2.*1.042
         j = INT(ekinetic)
         if (residue)THEN
             mult_10Be = mult_10be + weight
             ek_10be = ek_10be + ekinetic*weight
             ek2_10be = ek2_10be + ekinetic**2.*weight
         END IF
         IF (j <= 100) THEN
             ek(13,j) = ek(13,j) + weight
             IF (residue) ek2(13,j) = ek2(13,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 0)&
             ek3(13,j) = ek3(13,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 1)&
             ek4(13,j) = ek4(13,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(13,j) = ang(13,j)&
         + weight/SIN_D(frag(i)%vel%theta)


     ELSE IF (frag(i)%Z == 2 .AND. frag(i)%a == 6) THEN
         ekinetic = 0.5*6.*frag(i)%vel%v**2.*1.042
         j = INT(ekinetic)
         if (residue)THEN
             mult_6He = mult_6he + weight
             ek_6he = ek_6he + ekinetic*weight
             ek2_6he = ek2_6he + ekinetic**2.*weight
         END IF
         IF (j <= 100) THEN
             ek(11,j) = ek(11,j) + weight
             IF (residue) ek2(11,j) = ek2(11,j) + weight
             if (residue .AND. i == 1) THEN
                 ek4(11,j) = ek4(11,j) + weight
             ELSE if (residue .AND.&
              frag(i)%id /= frag(i-1)%id) THEN
                 ek4(11,j) = ek4(11,j) + weight
             ELSE
                 ek3(11,j) = ek3(11,j) + weight
             END IF
          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(11,j) = ang(11,j)&
         + weight/SIN_D(frag(i)%vel%theta)

     ELSE IF (frag(i)%Z == 4 .AND. frag(i)%a == 9) THEN
         ekinetic = 0.5*9.*frag(i)%vel%v**2.*1.042
         j = INT(ekinetic)
         if (residue)THEN
             mult_9Be = mult_9be + weight
             ek_9be = ek_9be + ekinetic*weight
             ek2_9be = ek2_9be + ekinetic**2.*weight
         END IF
         IF (j <= 100) THEN
             ek(12,j) = ek(12,j) + weight
             IF (residue) ek2(12,j) = ek2(12,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 0)&
             ek3(12,j) = ek3(12,j) + weight
             if (residue .AND. MODULO(frag(i)%id,2) == 1)&
             ek4(12,j) = ek4(12,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(12,j) = ang(12,j)&
         + weight/SIN_D(frag(i)%vel%theta)

     ELSE IF (frag(i)%Z == 3 .OR. frag(i)%Z == 4) THEN
       vzz = frag(i)%vel%v*COS_D(frag(i)%vel%theta) + vs
       vxx = frag(i)%vel%v*SIN_D(frag(i)%vel%theta)
       vv = (vxx**2.+vzz**2.)**0.5
      if (vv > 0) then
         theta_lab_Z = ACOS_D(vzz/vv)
      else
         theta_lab_Z = 0
         WRITE (UNIT=*, FMT=*) "li,Be"
      END IF

       ekinetic = 0.5*frag(i)%a*vv**2.*1.042
         j = INT(ekinetic/2)
         k = theta_lab_Z/20 + 1
         ii  = frag(i)%Z - 2
         IF (j <= 100) THEN
             IF (residue) ek_h(ii,j,k) = ek_h(ii,j,k) + weight
         END IF
      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)%a > 0 .AND. frag(i)%Z > 0) THEN
!            increment Charge distributions
     IZ = frag(I)%Z + 1
     A_dist(frag(i)%a) = A_dist(frag(i)%a) + weight
     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.

!    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
       ELSE IF (iZmax2 < frag(I)%Z) THEN
         iZmax2 = frag(i)%Z
       END IF
     END IF
   END IF
END DO





 IF (frag(n_frag)%Z >= z_ev) THEN
   j = theta_lab*2
   IF (j <= 50) THEN
       IF (theta_lab <= .05) theta_lab = .05
       ang_res(j,-1) = ang_res(j,-1) + weight/SIN_D(theta_lab)
       if (alpha <= 5) THEN
          ang_res(j,alpha) = ang_res(j,alpha) + weight/SIN_D(theta_lab)
          xsec_er(alpha) = xsec_er(alpha) + weight
       END IF
   END IF
 END IF


 if (residue .and. Mult_aa > 1) THEN
     DO i=1,mult_aa-1
         DO j=i+1,mult_aa
             delta_phi = phi_a(i) - phi_a(j)
             if (delta_phi < 0) delta_phi = delta_phi + 360
             ix = delta_phi/10
             phi_aa(ix) = phi_aa(ix) + weight
             ix =(360-delta_phi)/10
             phi_aa(ix) = phi_aa(ix) + weight
         END DO
     END DO
 END IF


!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
RETURN


END SUBROUTINE STATISTICS


SUBROUTINE OUTPUT ()
!
! FUNCTIONAL DESCRIPTION:
!
!    This subroutine writes out the contents of the histograms
!
!
!    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
INTEGER :: k,i,j
REAL (KIND=r4) :: dist,dist_2
!REAL (KIND=r4) :: da,sumation,mass_asy


REAL (KIND=r4) :: mean_e,sig_e

IF (N_residue > 0) THEN
write (UNIT=8,FMT=*) " <Zer> = ",sumzer/n_residue
write (UNIT=8,FMT=*) " <Ner> = ",sumner/n_residue
END IF

if (n_all > 0) WRITE (UNIT=8, FMT=*) "ER inside of PPAC=",n_ins/n_all


WRITE (UNIT=*, FMT=*) "RESIDUES=",N_RESIDUE/para%r_num*factor,&
"+-",N_RESIDUE**0.5/para%r_num*factor, "mb"
WRITE (UNIT=8, FMT=*)"RESIDUES=",N_RESIDUE/para%r_num*factor,&
"+-",N_RESIDUE**0.5/para%r_num*factor, "mb"

WRITE (*,*) "multiplicity,  error, average energy, error" 
if (mult_n >= 2) THEN
  mean_e = ek_n/mult_n
  sig_e = ek2_n/(mult_n-1)
  sig_e = sig_e - mult_n/(mult_n-1)*mean_e**2.
  sig_e = sig_e**0.5
  sig_e = sig_e/mult_n**0.5
  WRITE (UNIT=8, FMT=*) "neutron ",&
  Mult_n/n_residue,mult_n**0.5/n_residue,mean_e,sig_e
END IF

if (mult_p >= 2) THEN
  mean_e = ek_p/mult_p
  sig_e = ek2_p/(mult_p-1)
  sig_e = sig_e - mult_p/(mult_p-1)*mean_e**2.
  sig_e = sig_e**0.5
  sig_e = sig_e/mult_p**0.5
  WRITE (UNIT=8, FMT=*) "proton ",&
  Mult_p/n_residue,mult_p**0.5/n_residue,mean_e,sig_e
END IF

if (mult_d >= 2) THEN
  mean_e = ek_d/mult_d
  sig_e = ek2_d/(mult_d-1)
  sig_e = sig_e - mult_d/(mult_d-1)*mean_e**2.
  sig_e = sig_e**0.5
  sig_e = sig_e/mult_d**0.5
  WRITE (UNIT=8, FMT=*) "deuteron",&
  Mult_d/n_residue,mult_d**0.5/n_residue,mean_e,sig_e
END IF

if (mult_t >= 2) THEN
  mean_e = ek_t/mult_t
  sig_e = ek2_t/(mult_t-1)
  sig_e = sig_e - mult_t/(mult_t-1)*mean_e**2.
  sig_e = sig_e**0.5
  sig_e = sig_e/mult_t**0.5
  WRITE (UNIT=8, FMT=*) "triton",&
  Mult_t/n_residue,mult_t**0.5/n_residue,mean_e,sig_e
END IF

if (mult_a >= 2) THEN
  mean_e = ek_a/mult_a
  sig_e = ek2_a/(mult_a-1)
  sig_e = sig_e - mult_a/(mult_a-1)*mean_e**2.
  sig_e = sig_e**0.5
  sig_e = sig_e/mult_a**0.5
  WRITE (UNIT=8, FMT=*) "alpha",&
  Mult_a/n_residue,mult_a**0.5/n_residue,mean_e,sig_e
END IF

if (mult_3he >= 2) THEN
  mean_e = ek_3he/mult_3he
  sig_e = ek2_3he/(mult_3he-1)
  sig_e = sig_e - mult_3he/(mult_3he-1)*mean_e**2.
  sig_e = sig_e**0.5
  sig_e = sig_e/mult_3he**0.5
  WRITE (UNIT=8, FMT=*) "3he",&
  Mult_3he/n_residue,mult_3he**0.5/n_residue,mean_e,sig_e
END IF

if (mult_6he >= 2) THEN
  mean_e = ek_6he/mult_6he
  sig_e = ek2_6he/(mult_6he-1)
  sig_e = sig_e - mult_6he/(mult_6he-1)*mean_e**2.
  sig_e = sig_e**0.5
  sig_e = sig_e/mult_6he**0.5
  WRITE (UNIT=8, FMT=*) "6He=",&
  Mult_6he/n_residue,mult_6he**0.5/n_residue,mean_e,sig_e
END IF

if (mult_6li >= 2) THEN
  mean_e = ek_6li/mult_6li
  sig_e = ek2_6li/(mult_6li-1)
  sig_e = sig_e - mult_6li/(mult_6li-1)*mean_e**2.
 if (sig_e < 0) sig_e = 0
  sig_e = sig_e**0.5
  sig_e = sig_e/mult_6li**0.5
  WRITE (UNIT=8, FMT=*) "6Li=",&
  Mult_6li/n_residue,mult_6li**0.5/n_residue,mean_e,sig_e
END IF

if (mult_7li >= 2) THEN
  mean_e = ek_7li/mult_7li
  sig_e = ek2_7li/(mult_7li-1)
  sig_e = sig_e - mult_7li/(mult_7li-1)*mean_e**2.
 if (sig_e < 0) sig_e = 0
  sig_e = sig_e**0.5
  sig_e = sig_e/mult_7li**0.5
  WRITE (UNIT=8, FMT=*) "7Li",&
  Mult_7li/n_residue,mult_7li**0.5/n_residue,mean_e,sig_e
END IF

if (mult_8li >= 2) THEN
  mean_e = ek_8li/mult_8li
  sig_e = ek2_8li/(mult_8li-1)
  sig_e = sig_e - mult_8li/(mult_8li-1)*mean_e**2.
 if (sig_e < 0) sig_e = 0
  sig_e = sig_e**0.5
  sig_e = sig_e/mult_8li**0.5
  WRITE (UNIT=8, FMT=*) "8Li",&
  Mult_8li/n_residue,mult_8li**0.5/n_residue,mean_e,sig_e
END IF

if (mult_7be >= 2) THEN
  mean_e = ek_7be/mult_7be
  sig_e = ek2_7be/(mult_7be-1)
  sig_e = sig_e - mult_7be/(mult_7be-1)*mean_e**2.
 if (sig_e < 0) sig_e = 0
  sig_e = sig_e**0.5
  sig_e = sig_e/mult_7be**0.5
  WRITE (UNIT=8, FMT=*) "7Be",&
  Mult_7be/n_residue,mult_7be**0.5/n_residue,mean_e,sig_e
END IF

if (mult_8be >= 2) THEN
  mean_e = ek_8be/mult_8be
  sig_e = ek2_8be/(mult_8be-1)
  sig_e = sig_e - mult_8be/(mult_8be-1)*mean_e**2.
 if (sig_e < 0) sig_e = 0
  sig_e = sig_e**0.5
  sig_e = sig_e/mult_8be**0.5
  WRITE (UNIT=8, FMT=*) "8Be",&
  Mult_8be/n_residue,mult_8be**0.5/n_residue,mean_e,sig_e
END IF

if (mult_9be >= 2) THEN
  mean_e = ek_9be/mult_9be
  sig_e = ek2_9be/(mult_9be-1)
  sig_e = sig_e - mult_9be/(mult_9be-1)*mean_e**2.
 if (sig_e < 0) sig_e = 0
  sig_e = sig_e**0.5
  sig_e = sig_e/mult_9be**0.5
  WRITE (UNIT=8, FMT=*) "9Be",&
  Mult_9be/n_residue,mult_9be**0.5/n_residue,mean_e,sig_e
END IF

if (mult_10be >= 2) THEN
  mean_e = ek_10be/mult_10be
  sig_e = ek2_10be/(mult_10be-1)
  sig_e = sig_e - mult_10be/(mult_10be-1)*mean_e**2.
 if (sig_e < 0) sig_e = 0
  sig_e = sig_e**0.5
  sig_e = sig_e/mult_10be**0.5
  WRITE (UNIT=8, FMT=*) "10Be",&
  Mult_10be/n_residue,mult_10be**0.5/n_residue,mean_e,sig_e
END IF

If (sum_gam_energy > 0. )THEN
  write (UNIT=8, FMT=*) "total gamma energy = ", sum_gam_energy/n_residue
  write (UNIT=8, FMT=*) "missing energy(GAM+ER_EK) = ", sum_missing/n_residue
END IF

WRITE (UNIT=6, FMT="(TR2,A17,TR1,F9.2,A2)")&
"Fusion Xsection=",fusion*factor/para%r_num,"mb"
WRITE (UNIT=8, FMT="(TR2,A17,TR1,F9.2,A2)")&
"Fusion Xsection=",fusion*factor/para%r_num,"mb"

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/para%r_num
   dist_2 = Z_dist_2(i)**0.5*factor/para%r_num
   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 = 200
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/para%r_num
   WRITE (UNIT=8, FMT="(TR2,I3,TR2,F19.13,TR2,F19.13,TR2,f8.0)")&
   I, DIST
   END IF
END DO


!      write coincidence results
!       WRITE (UNIT=8, FMT="(/TR3,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", "<Zloss>", "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) = 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", "<A>", "sig(A)", "<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)<=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 == 1) THEN
!         WRITE (UNIT=8, FMT="(/TR2,a,/TR3,a1,TR7,a5,TR6,a8)")&
!         "Average c-m velocity of each Z-species",&
!         "Z", "<Vcm>", "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<Z<21)",&
!           "for events with at least one IMF=", akp3
!          END IF

!       WRITE (UNIT=8, FMT="(TR2,a35)")&
!       "pre and post fission multiplicities"
!       WRITE (UNIT=8, FMT="(TR2,a4,TR13,a3,TR14,a15,TR8,a15)")&
!       "mass","pre","post-heavy-frag","post-light-frag"
!       WRITE (UNIT=8, FMT="(TR2,a3,3(TR2,3(TR4,a3)))")&
!       "asy","n","p","alp","n","p","alp","n","p","alp"
!       DO i=0,9
!           IF (total(i) > 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)
!               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)
!             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"

WRITE (UNIT=8, FMT=*)
WRITE (unit=8, FMT="(TR2,A)")&
"c.m. kinetic energy (MeV), dsigma/de (mb/MeV) for n,p,a"
IF (para%i_angle) THEN
    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

if (n_residue > 0) THEN
   WRITE (UNIT=8, FMT=*)
   WRITE (unit=8, FMT="(TR2,A)")&
   "c.m. kinetic energy (MeV), differential multiplicity dm/de/domega (1/MeV/sr) 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)")&
   "c.m. kinetic energy (MeV), dm/de/domega (1/MeV/sr)  n,p,a from sequential decay"

    DO j=0,100
        WRITE (UNIT=8, FMT="(TR2,F5.1,3(TR2,E10.5))")&
        REAL(j)+.5,ek3(1,j)/n_residue/12.56,&
        ek3(2,j)/n_residue/12.56,ek3(3,j)/n_residue/12.56
     END DO

   WRITE (UNIT=8, FMT=*)
   WRITE (unit=8, FMT="(TR2,A)")&
   "c.m. kinetic energy (MeV), dm/de/domega  n,p,a from non-sequential decay"

    DO j=0,100
        WRITE (UNIT=8, FMT="(TR2,F5.1,3(TR2,E10.5))")&
        REAL(j)+.5,ek4(1,j)/n_residue/12.56,&
        ek4(2,j)/n_residue/12.56,ek4(3,j)/n_residue/12.56
    END DO


   WRITE (UNIT=8, FMT=*)
   WRITE (unit=8, FMT="(TR2,A)")&
   "c.m. angular distribtion relative to spin axis n,p,a"
   WRITE (UNIT=8, FMT=*) "  angle (deg), dm/domega (1/sr)"
    DO j=0,17
        WRITE (UNIT=8, FMT="(TR2,I3,3(TR2,E10.5))")&
        j*10+5,ang(1,j)/10*180/3.14159/n_residue/6.2832&
        ,ang(2,j)/10*180/3.14159/N_residue/6.2832&
        ,ang(3,j)/10*180/3.14159/n_residue/6.2832
     END DO
END IF

   WRITE (UNIT=8, FMT=*)
   WRITE (unit=8, FMT="(TR2,A)")&
   "c.m. kinetic energy, d_sigma/de (mb/MeV) 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

if (n_residue > 0) THEN
      WRITE (UNIT=8, FMT=*)
      WRITE (unit=8, FMT="(TR2,A)")&
      "c.m. kinetic energy, dm/de/domega 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)")&
      "c.m. angular distribtion dm/domega (1/sr) 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/n_residue/6.2832&
           ,ang(5,j)/10*180/3.14159/N_residue/6.2832&
           ,ang(6,j)/10*180/3.14159/n_residue/6.2832
       END DO



      WRITE (UNIT=8, FMT=*)
      WRITE (unit=8, FMT="(TR2,A)")&
      "c.m. kinetic energy, dm/de/domega (1/MeV/sr) for 6li,7li,8li from residue"

       DO j=0,100
           WRITE (UNIT=8, FMT="(TR2,F5.1,3(TR2,E10.5))")&
           REAL(j)+.5,ek2(7,j)/n_residue/12.56,&
           ek2(8,j)/n_residue/12.56,&
           ek2(9,j)/n_residue/12.56
       END DO


      WRITE (UNIT=8, FMT=*)
      WRITE (unit=8, FMT="(TR2,A)")&
      "c.m. angular distribtion dm/domega (1/sr) relative to spin axis 6li,7li,li"
       DO j=0,17
           WRITE (UNIT=8, FMT="(TR2,I3,3(TR2,E10.5))")&
           j*10+5,ang(7,j)/10*180/3.14159/n_residue/6.2832&
           ,ang(8,j)/10*180/3.14159/N_residue/6.2832&
           ,ang(9,j)/10*180/3.14159/n_residue/6.2832
       END DO


      WRITE (UNIT=8, FMT=*)
      WRITE (unit=8, FMT="(TR2,A)")&
      "c.m. kinetic energy, dm/de/domega (1/MeV/sr) for 6He,7Be,9be from residue"

       DO j=0,100
           WRITE (UNIT=8, FMT="(TR2,F5.1,3(TR2,E10.5))")&
           REAL(j)+.5,ek2(11,j)/n_residue/12.56,&
           ek2(10,j)/n_residue/12.56,&
           ek2(12,j)/n_residue/12.56
       END DO


      WRITE (UNIT=8, FMT=*)
      WRITE (unit=8, FMT="(TR2,A)")&
      "c.m. angular distribtion dm/domega (1/sr) relative to spin axis 6He,7Be,9Be"
       DO j=0,17
           WRITE (UNIT=8, FMT="(TR2,I3,3(TR2,E10.5))")&
           j*10+5,ang(11,j)/10*180/3.14159/n_residue/6.2832&
           ,ang(10,j)/10*180/3.14159/N_residue/6.2832&
           ,ang(12,j)/10*180/3.14159/N_residue/6.2832
       END DO


      WRITE (UNIT=8, FMT=*)
      WRITE (unit=8, FMT="(TR2,A)")&
      "c.m. kinetic energy, dm/de/domega (1/MeV/sr) for 10be,8Be,8be primary from residue"

       DO j=0,100
           WRITE (UNIT=8, FMT="(TR2,F5.1,3(TR2,E10.5))")&
           REAL(j)+.5,ek2(13,j)/n_residue/12.56,&
           ek2(14,j)/n_residue/12.56,ek4(14,j)/n_residue/12.56
       END DO


      WRITE (UNIT=8, FMT=*)
      WRITE (unit=8, FMT="(TR2,A)")&
      "c.m. angular distribtion dm/domega relative to spin axis 10Be,8Be"
       DO j=0,17
           WRITE (UNIT=8, FMT="(TR2,I3,3(TR2,E10.5))")&
           j*10+5,ang(13,j)/10*180/3.14159/n_residue/6.2832&
           ,ang(13,j)/10*180/3.14159/N_residue/6.2832
       END DO
   END IF

   WRITE (unit=8, FMT=*) " residue angular distribution gated on the the alpha multiplicity"

   WRITE (UNIT=8, FMT=*)
   write (UNIT=8, FMT="(TR2,A11,6(TR2,E9.4))")&
   "xsec (ER)=",xsec_er(0)*factor/para%r_num,&
   xsec_er(1)*factor/para%r_num,xsec_er(2)*factor/para%r_num,&
   xsec_er(3)*factor/para%r_num,xsec_er(4)*factor/para%r_num,&
   xsec_er(5)*factor/para%r_num
   WRITE (unit=8, FMT="(TR2,A)")&
   " residue angular distribtion dsigma/domega in mb/sr, mult= alpha-particle multiplicity"
   WRITE (UNIT=8, FMT=*) " angle (deg), mult=0, mult=2, mult=3, mult=4,  mult=5 "
    DO j=0,50
        WRITE (UNIT=8, FMT="(TR2,F5.2,6(TR2,E9.4))")&
        (REAL(j)+.5)/2&
        ,ang_res(j,0)/6.283/3.14159*180.*factor/para%r_num*2&
        ,ang_res(j,1)/6.283/3.14159*180.*factor/para%r_num*2&
        ,ang_res(j,2)/6.283/3.14159*180.*factor/para%r_num*2&
        ,ang_res(j,3)/6.283/3.14159*180.*factor/para%r_num*2&
        ,ang_res(j,4)/6.283/3.14159*180.*factor/para%r_num*2&
        ,ang_res(j,5)/6.283/3.14159*180.*factor/para%r_num*2
        END DO



    WRITE (UNIT=8, FMT=*)
    WRITE (UNIT=8, FMT=*) "  all residues"
    WRITE (UNIT=8, FMT=*) "  angle (deg), dsigma/domega (mb/sr)"
    DO j=0,50
        WRITE (UNIT=8, FMT="(TR2,F5.2,6(TR2,E9.4))")&
        (REAL(j)+.5)/2&
        ,ang_res(j,-1)/6.283/3.14159*180.*factor/para%r_num*2
    END DO


!            DO j=2,7
!               WRITE (UNIT=89, FMT=*) ang_res_step(j)*factor
!           END DO


     if (n_residue > 0) THEN
       WRITE (UNIT=8, FMT=*)
       WRITE (unit=8, FMT="(TR2,A)")&
       "residue velocity distribtion gated on selected angular range theta_min-theta_max"


        WRITE (UNIT=8, FMT=*) " lab velocity (cm/ns), dsigma/dv (mb/cm/ns)"
        WHERE (v_residue /= 0.)
           av_residue = av_residue/v_residue
        END WHERE
        DO j=0,120
             WRITE (UNIT=8, FMT="(TR2,F5.3,TR2,E10.5,tr2,F6.2)")&
             (REAL(j)+.5)/40,v_residue(j)/n_residue,av_residue(j)
        END DO

       WRITE (UNIT=8, FMT=*)
       WRITE (unit=8, FMT="(TR2,A)")&
       "angle gated residue mass distribtion"
       WRITE (UNIT=8, FMT=*) "nucleon number (Mass) , sigma (mb)"
        DO j=50,220
            WRITE (UNIT=8, FMT="(TR2,I3,TR2,E10.5,TR2,E10.5)")&
            j,a_residue(j)*factor/para%r_num,&
            a_residue(j)**0.5*factor/para%r_num
        END DO
      END IF


       WRITE (UNIT=8, FMT=*)
       WRITE (unit=8, FMT="(TR2,A)")&
       "a-a phi-phi correlation"
    DO j=0,35
        WRITE (UNIT=8, FMT="(TR2,F6.2,TR2,E10.5)")&
        (REAL(j)+.5)*10,phi_aa(j)
        END DO

    END IF



!       OPEN (
!       1 UNIT=55,
!       1 FILE = "li.out",
!       1 STATUS = "NEW")
!       DO k=1,9
!           DO j=1,100
!               WRITE (UNIT=55, FMT=*) 2*j+1,ek_h(1,j,k)
!           END DO
!       END DO
!
!       OPEN (
!       1 UNIT=56,
!       1 FILE = "be.out",
!       1 STATUS = "NEW")
!       DO k=1,9
!           DO j=1,100
!               WRITE (UNIT=56, FMT=*) 2*j+1,ek_h(2,j,k)
!           END DO
!       END DO






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
!
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

END PROGRAM sort_er