PROGRAM LOOK
!
! PROGRAM DESCRIPTION:
!
!      prints out decay chains of gemimi events
!      to continue reading out events after fortran pause type
!      "continue" or just  'c'.
!
! AUTHORS:
!
!      Bob Charity (RJC)
!
! CREATION DATE:        1988
!
!
!          C H A N G E   L O G
!
!      Date     | Name  | Description
! ----------------+-------+-----------------------------------------------------
! [change_entry]
!
USE gem_kind
USE trig
USE sortsub
IMPLICIT NONE


INTEGER :: L
INTEGER :: I_read_error,N_tot,M_frag,I,jj,zzz,j,id
INTEGER :: unit_in,unit_out,events,ie
INTEGER, PARAMETER :: n_hits=150
INTEGER, DIMENSION(n_hits) :: hits
REAL (KIND=r4) :: sumaa2,sumnz,sumnz2,sumzat,sumzat2,zat
REAL (KIND=r4) :: sump,sumd,sumt,sumhe3,suma,sumhe6,sumli,summ
REAL (KIND=r4), DIMENSION(0:150) :: sums,sume
REAL (KIND=r4) :: jmax,dif,sumz,sumaa,sumbe7,sumbe9,sumbe10
REAL (KIND=r4) :: sumbe11,sumn,sum_binary,sumbe
REAL (KIND=r4) :: vx1,vx2,vy1,vy2,vz1,vz2,e_rel,delta_v
REAL (KIND=r4) :: sumbe8,sumbe81,sumli6,sumli7,sumli8,sumli9
REAL (KIND=r4) :: sumli7_2,sumli6_1,sumhe5,sumhe7,sumhe8
REAL (KIND=r4), DIMENSION(0:50) :: spectra_he5,spectra_li6
REAL (KIND=r4) :: tot_z,tot_a,sumbe9_2,sumbe7_1,sumbe8_p
REAL (KIND=r4) :: sumhe5_p,sumli5,sumli5_p
REAL (KIND=r4) :: sum_gam_energy  
LOGICAL :: primary
LOGICAL, SAVE :: first=.true.


sums = 0.0
sume = 0.0
sum_binary = 0.0
sumaa = 0.0
summ = 0.0
sumz = 0.0
sumaa2 = 0.0
sumnz =  0.0
sumnz2 = 0.0
sumzat = 0.0
sumzat2 = 0.0
spectra_he5 = 0.0
spectra_li6 = 0.0
events = 0
sumn = 0.0
sump = 0.0
sumd = 0.0
sumt = 0.0
sumhe3 = 0.0
suma = 0.0
sumhe5 = 0.0
sumhe6 = 0.0
sumhe7 = 0.0
sumhe8 = 0.0
sumli = 0.0
sumli6 = 0.0
sumli7 = 0.0
sumli8 =0.0
sumli9 = 0.0
sumli6_1 = 0.0
sumli7_2 = 0.0
sumbe = 0.0
sumbe7 = 0.0
sumbe7_1 = 0.0
sumbe8 = 0.0
sumbe81 = 0.0
sumbe9 = 0.0
sumbe9_2 = 0.0
sumbe10 = 0.0
sumbe11 = 0.0
sumbe8_p =0.0
sumhe5_p = 0.0
sumli5 = 0.0
sumli5_p = 0.0
sum_gam_energy = 0.0

N_TOT = 0
M_frag = 0


DO

    unit_in = 2
    unit_out = -1



    call start_sort(unit_in,unit_out,i_read_error)
    if (I_read_error /= 0) EXIT
    if (first) THEN
        WRITE (UNIT=*, FMT=*) "z_min,jmax,dif"
        READ (UNIT=5, FMT=*) zzz,jmax,dif
        first = .FALSE.
    END IF 

    DO 
       call get_event(l,i_read_error)
       if (I_read_error /= 0)EXIT
       events = events + 1

       IF (dif == 0) THEN
           IF (l > jmax) weight = 0.0
       ELSE
           weight = weight/(1.0+exp((REAL(l)-jmax)/dif))
       END IF

       sums(l) = sums(l) + weight
       IF (frag(n_frag)%Z >= zzz) THEN
            sume(l) = sume(l) + weight
            sum_gam_energy = sum_gam_energy + weight*frag(n_frag)%exx
            sum_binary = sum_binary + REAL(N_binary)*weight
            summ = summ  + weight
            sumz = sumz + weight*frag(n_frag)%z
            sumaa = sumaa + weight*frag(n_frag)%a
            sumaa2 = sumaa2  + weight*frag(n_frag)%a**2
            sumnz = sumnz + weight*(frag(n_frag)%a-2*frag(n_frag)%z)
            sumnz2 = sumnz2 + weight*(frag(n_frag)%a-2*frag(n_frag)%z)**2.0
            zat = 0.479*frag(n_frag)%A - 2.04e-4*REAL(frag(n_frag)%A)**2.0
            sumzat = sumzat + weight*(REAL(frag(n_frag)%z)-Zat)
            sumzat2 = sumzat2 + weight*(REAL(frag(n_frag)%z)-Zat)**2.0
            tot_z = frag(n_frag)%z
            tot_a = frag(n_frag)%a

            hits = 0

            DO i=1,n_frag-1
                id = frag(i)%id
		if (id > n_hits) THEN
                   WRITE (*,*) "increase n_hits in prob.f90"
		   STOP
		END IF 
                hits(id) = hits(id) + 1


                if (frag(i)%z < 0 .AND. frag(i)%A <= 2)THEN
!                  gamma ray
                   sum_gam_energy = sum_gam_energy + weight*frag(i)%exx
		END IF 
                if (frag(i)%z >= 0 .AND. frag(i)%A > 0)THEN
                    tot_z = tot_z + frag(i)%z
                    tot_a = tot_a + frag(i)%a
                END IF 

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

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

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

                    if (frag(i)%A == 4 .AND. frag(i-1)%a == 4) THEN

                    e_rel = delta_v/0.9784**2.0*2.0/2.0
                    if (e_rel <= 0.15) THEN
                        sumbe8 = sumbe8 + weight
                        primary = .TRUE.
                        DO jj=1,i-2
                            if (frag(jj)%id == frag(i)%id) THEN
                                primary = .FALSE.
                                EXIT
                             END IF 
                        END DO 
  
                        if (primary) sumbe8_p = sumbe8_p + weight
 
                    ELSE
                        sumbe81 = sumbe81 + weight
                    END IF

                   ELSE IF (frag(i)%A == 4&
                   .AND. frag(i-1)%a == 3 .AND. frag(i-1)%Z == 1) THEN
                      sumli7_2 = sumli7_2 + weight

                   ELSE IF (frag(i)%A == 4&
                   .AND. frag(i-1)%a == 3 .AND. frag(i-1)%Z == 2) THEN
                      sumbe7_1 = sumbe7_1 + weight

                   ELSE IF (frag(i)%A == 4&
                   .AND. frag(i-1)%a == 2 .AND. frag(i-1)%Z == 1) THEN
                      e_rel = delta_v/0.9784**2.0*8.0/6.0/2.0
                      if (e_rel <= 0.8)sumli6_1 = sumli6_1 + weight
    !                           sumli6_1 = sumli6_1 + weight
                      ie = INT(e_rel*5)
                      spectra_li6(ie) = spectra_li6(ie) + weight

                  ELSE IF (frag(i)%A == 4 .AND. frag(i)%Z == 2&
                  .AND. frag(i-1)%a == 1 .AND. frag(i-1)%Z == 1)THEN
                     if (hits(frag(i)%id) == 2) THEN

                       e_rel = delta_v/0.9784**2.0*4.0/5.0/2.0
                       ie = INT(e_rel*5.0)
                       sumli5 = sumli5 + weight
                      END IF   


                  ELSE IF (frag(i)%A == 4 .AND. frag(i)%Z == 2&
                 .AND. frag(i-1)%a == 1 .AND. frag(i-1)%Z == 0)THEN
                      if (hits(frag(i)%id) == 2) THEN

                       e_rel = delta_v/0.9784**2.0*4.0/5.0/2.0
                       ie = INT(e_rel*5.0)
                       sumhe5 = sumhe5 + weight

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

                       if (ie <= 50)spectra_he5(ie)= spectra_he5(ie) + weight
                       else if (hits(frag(i)%id) == 3) THEN
                       j = i-2
                       DO 
                          if (frag(j)%id == frag(i)%id) EXIT
                              j = j - 1
                          END DO
                          if (frag(j)%A == 4) THEN
                              sumbe9_2 = sumbe9_2  + weight
                          END IF 
                       END IF 

                  ELSE IF (frag(i)%A == 6 .AND. frag(i)%Z == 2&
                  .AND. frag(i-1)%a == 1 .AND. frag(i-1)%Z == 0) THEN
                      sumhe7 = sumhe7 + weight
                  end if
                END IF 
                END IF 


                IF (frag(i)%Z == 0 .AND. frag(i)%A == 1) THEN
                    sumn = sumn + weight
                ELSE IF (frag(i)%Z == 1 .AND. frag(i)%A == 1) THEN
                    sump = sump + weight
                ELSE IF (frag(i)%Z == 1 .AND. frag(i)%A == 2) THEN
                    sumd = sumd + weight
                ELSE IF (frag(i)%Z == 1 .AND. frag(i)%A == 3) THEN
                    sumt = sumt + weight
                ELSE IF (frag(i)%Z == 2 .AND. frag(i)%A == 3) THEN
                    sumhe3 = sumhe3+ weight
                ELSE IF (frag(i)%Z == 2 .AND. frag(i)%A == 4) THEN
                    suma = suma + weight
                ELSE IF (frag(i)%Z == 2 .AND. frag(i)%A == 6) THEN
                    sumhe6 = sumhe6 + weight
                ELSE IF (frag(i)%Z == 2 .AND. frag(i)%A == 8) THEN
                    sumhe8 = sumhe8 + weight
                ELSE IF (frag(I)%Z == 3) THEN
                    sumli = sumli + weight
                    if (frag(i)%A == 6) THEN
                        sumli6 = sumli6 + weight
                    ELSE if (frag(i)%A == 7) THEN
                        sumli7 = sumli7 + weight
                    ELSE if (frag(i)%A == 8) THEN
                        sumli8 = sumli8 + weight
                    ELSE if (frag(i)%A == 9) THEN
                        sumli9 = sumli9 + weight
                    END IF 

               ELSE IF (frag(I)%Z == 4) THEN
                    sumbe= sumbe + weight
                    if (frag(i)%A == 7) THEN
                        sumbe7 = sumbe7 + weight
                    ELSE if (frag(i)%A == 9) THEN
                        sumbe9 = sumbe9 + weight
                    ELSE if (frag(i)%A == 10) THEN
                        sumbe10 = sumbe10 + weight
                    ELSE if (frag(i)%A == 11) THEN
                        sumbe11 = sumbe11 + weight
                    END IF 
               END IF 
            END DO
            if (tot_z /= para%zcn .AND. tot_a /= para%Acn)THEN
              WRITE (UNIT=*, FMT=*) tot_z,tot_a
            END IF 
        END IF
    END DO
END DO
CLOSE (UNIT=2)
DO i=0,150
    IF (sums(i) > 0) THEN
        WRITE (UNIT=54, FMT=*) i,sume(i)/sums(i)
        WRITE (UNIT=55, FMT=*) i,sume(i)
    END IF
END DO

WRITE (UNIT=*, FMT=*) "events=",events
write (*,*) 'residue events=',summ
WRITE (UNIT=*, FMT=*) "neutron=",sumn/summ,sumn**0.5/summ,sumn
WRITE (UNIT=*, FMT=*) "proton =",sump/summ,sump**0.5/summ,sump
WRITE (UNIT=*, FMT=*) "deutron=",sumd/summ,sumd**0.5/summ,sumd
WRITE (UNIT=*, FMT=*) "triton =",sumt/summ,sumt**0.5/summ,sumt
WRITE (UNIT=*, FMT=*) "3He    =",sumhe3/summ,sumhe3**0.5/summ,sumhe3
WRITE (UNIT=*, FMT=*) "alpha  =",suma/summ,suma**0.5/summ,suma
WRITE (UNIT=*, FMT=*) "5He    =",sumhe5/summ,sumhe5**0.5/summ,sumhe5
WRITE (UNIT=*, FMT=*) "5he gs pramary=",sumhe5_p/summ,sumhe5_p**0.5/summ,sumhe5_p
if (sumhe6 > 0) &
WRITE (UNIT=*, FMT=*) "6He    =",sumhe6/summ,sumhe6**0.5/summ,sumhe6
if (sumhe7 > 0) &
WRITE (UNIT=*, FMT=*) "7He    =",sumhe7/summ,sumhe7**0.5/summ,sumhe7
if (sumhe8 > 0) &
WRITE (UNIT=*, FMT=*) "8He    =",sumhe8/summ,sumhe8**0.5/summ,sumhe8

if (sumli > 0) &
WRITE (UNIT=*, FMT=*) "Li     =",sumli/summ**0.5/summ,sumli
if (sumli5 > 0) &
WRITE (UNIT=*, FMT=*) "5Li     =",sumli5/summ,sumli5**0.5/summ,sumli5
if (sumli6 > 0) &
WRITE (UNIT=*, FMT=*) "6Li     =",sumli6/summ,sumli6**0.5/summ,sumli6
if (sumli6_1 > 0) &
WRITE (UNIT=*, FMT=*) "6Li(1st)=",sumli6_1/summ,sumli6_1**0.5/summ,sumli6_1
if (sumli7 > 0) &
WRITE (UNIT=*, FMT=*) "7Li     =",sumli7/summ,sumli7**0.5/summ,sumli7
if (sumli7_2 > 0) &
WRITE (UNIT=*, FMT=*) "7Li(2nd)=",sumli7_2/summ,sumli7_2**0.5/summ,sumli7_2
if (sumli8 > 0) &
WRITE (UNIT=*, FMT=*) "8Li     =",sumli8/summ,sumli8**0.5/summ,sumli8
if (sumli9 > 0) &
WRITE (UNIT=*, FMT=*) "9Li     =",sumli9/summ,sumli9**0.5/summ,sumli9

if (sumbe > 0) &
WRITE (UNIT=*, FMT=*) "BE     =",sumbe/summ**0.5/summ,sumbe
if (sumbe8 > 0) &
WRITE (UNIT=*, FMT=*) "8BEg.s =",sumbe8/summ,sumbe8**0.5/summ,sumbe8
if (sumbe8_p > 0) &
WRITE (UNIT=*, FMT=*) "8BE gs pramary=",sumbe8_p/summ,sumbe8_p**0.5/summ,sumbe8_p
if (sumbe81 > 0) &
WRITE (UNIT=*, FMT=*) "8BE(1st) =",sumbe81/summ,sumbe81**0.5/summ,sumbe81
if (sumbe7 > 0) &
WRITE (UNIT=*, FMT=*) "7BE     =",sumbe7/summ,sumbe7**0.5/summ,sumbe7
if (sumbe7_1 > 0) &
WRITE (UNIT=*, FMT=*) "7BE(1st)=",sumbe7_1/summ,sumbe7_1**0.5/summ,sumbe7_1
if (sumbe9 > 0) &
WRITE (UNIT=*, FMT=*) "9BE     =",sumbe9/summ,sumbe9**0.5/summ,sumbe9
if (sumbe9_2 > 0) &
WRITE (UNIT=*, FMT=*) "9BE(2nd)=",sumbe9_2/summ,sumbe9_2**0.5/summ,sumbe9_2
if (sumbe10 > 0) &
WRITE (UNIT=*, FMT=*) "10BE     =",sumbe10/summ,sumbe10**0.5/summ,sumbe10
if (sumbe11 > 0) &
WRITE (UNIT=*, FMT=*) "11BE     =",sumbe11/summ,sumbe11**0.5/summ,sumbe11


WRITE (UNIT=*, FMT=*) "primary imf=",sum_binary/summ
WRITE (UNIT=*, FMT=*) "SECONDARY IMF=",(SUMHE6+SUMLI+SUMBE)/SUMm
IF (sum_binary > 0)&
WRITE (UNIT=*, FMT=*) "SURVIVAL PROB=",(SUMHE6+SUMLI+SUMBE)/SUM_BINARY
WRITE (UNIT=*, FMT=*) "<Zer>=",sumz/summ
WRITE (UNIT=*, FMT=*) "<Aer>=",sumAa/summ
WRITE (UNIT=*, FMT=*) "<Ner>=",(sumAa-sumz)/summ
WRITE (UNIT=*, FMT=*) "<N-Z>=",sumnz/summ
sumaa2 = sumaa2/(summ-1.) - summ/(summ-1.)*(sumaa/summ)**2.0
if (sumaa2 <= 0.) THEN
     sumaa2 = 0.
ELSE
     sumaa2= sumaa2**0.5
END IF
WRITE (UNIT=*, FMT=*) "sigma A=",sumaa2
sumnz2 = (sumnz2/(summ-1.) - summ/(summ-1.)*(sumnz/summ)**2.0)**0.5
WRITE (UNIT=*, FMT=*) "sigma n-z=",sumnz2
WRITE (UNIT=*, FMT=*) "<z-Z_at>=",sumzat/summ
sumzat2 = (sumzat2/(summ-1.) - summ/(summ-1.)*(sumzat/summ)**2.0)**0.5
WRITE (UNIT=*, FMT=*) "sigma (z-z_at)=",sumzat2
! DO i=0,50
!     WRITE (UNIT=99, FMT=*) (REAL(i)+0.5)/5.0,spectra_he5(i)
! END DO
! DO i=0,50
!     WRITE (UNIT=98, FMT=*) (REAL(i)+0.5)/5.0,spectra_li6(i)
! END DO

write (UNIT=*,FMT=*) "average gam energy=", sum_gam_energy/summ

STOP
END PROGRAM LOOK