* Program to calculate differenced data for forward scattering
      CHARACTER*80 runid
      character runo*4,run1*5,fout*80
      REAL T(10000),CT(10000),err(10000),SUM1(9),SUM2(9),SUM3(9)
      REAL COUT(10000),CTHIN(10000),CTHICK(10000),NFR(9),sum4(9)
      REAL ERROUT(10000),ERRTHIN(10000),ERRTHICK(10000)
      REAL OUT(2000),THIN(2000),THICK(2000)
      REAL MOUT(10000),MTHIN(10000),MTHICK(10000)
      INTEGER IARRAY(100),ERRCODE,LOUT,NCHAN,NDET,IFSN,NOS
      INTEGER rpb1(32),ICT(10000),ICTA(9,10000),MICTA(9,10000),NI(9)
      real xin(10000),yin(10000),ein(10000) 
      real xout(10000),yout(10000),eout(10000) 

      write(6,*) ' First and last run numbers?'
      read(5,*) IRUN1,IRUN2 
      write(6,*) ' First and last spectrum numbers?'
      READ(5,*) ISPECT1,ISPECT2
*      WRITE(6,*) ' File name for output?'
*      read(5,'(a)') fout
      fout='temp.dat'
      open(unit=4,file=fout,status='new')
      
      DO 1000 ISPECT=ISPECT1,ISPECT2
      WRITE(6,*) ' spectrum: ',ISPECT
      IOPT=1
      IF(IOPT.EQ.1) THEN
       WRITE(6,*) ' THIN DIFFERENCE'
      ELSE IF(IOPT.EQ.2) THEN
       WRITE(6,*) ' DOUBLE DIFFERENCE'
       BETA=ST_VAR_IN(5)
       WRITE(6,*) ' BETA=',BETA
      ELSE IF(IOPT.EQ.3) THEN
       WRITE(6,*) ' THICK DIFFERENCE'
      END IF

      DO I=1,9
      DO J=1,10000
       ICTA(I,J)=0
       MICTA(I,J)=0
      END DO
      END DO 
      NGOOD_FRAMES=0    

      DO 101 IRUN=IRUN1,IRUN2
* Define name of data file.
      if(irun.lt.10000) then
       write(runo(1:4),'(I4.4)')irun
       runid='evs_data:evs0'//runo//'.raw'
      else
       write(run1(1:5),'(I5.5)')irun
       runid='evs_data:evs'//run1//'.raw'
      end if
*      write(6,5) runid
    5 Format(' ',' data file is;',50a)

* Get relevant run information.
      CALL GETPARI(runid,'RPB',rpb1,32,length_out,errcode)
      ngood_frames =ngood_frames+rpb1(10)
      CALL GETPARI(runid,'NPER',NPD,100,LOUT,ERRCODE)
      CALL GETPARI(RUNID,'NSP1',IARRAY,100,LOUT,ERRCODE)
      NDET=IARRAY(1)

      CALL GETPARR(RUNID,'TCB1',T,10000,LOUT,ERRCODE)
      CALL GETPARI(RUNID,'NTC1',IARRAY,100,LOUT,ERRCODE)
      NCHAN=IARRAY(1)
      NMON=MON(RUNID)  
      DT=T(3)-T(2)
      DO I=1,NCHAN
       T(I)=T(I)-DT
      END DO

      DO IP=1,NPD
       IS=ISPECT+(IP-1)*(NDET+1)
       CALL GETDAT(RUNID,IS,1,ICT,10000,ERRCODE)
       DO IT=2,NCHAN
        ICTA(IP,IT)=ICTA(IP,IT)+ICT(IT)
       END DO
       IS=NMON+(IP-1)*(NDET+1)
       CALL GETDAT(RUNID,IS,1,ICT,10000,ERRCODE)
       DO IT=2,NCHAN
        MICTA(IP,IT)=MICTA(IP,IT)+ICT(IT)
       END DO
      END DO

  101 CONTINUE

      write(6,*) ' Total number of good frames=',ngood_frames
      WRITE(6,*) ' NUMBER OF PERIODS=',NPD
      WRITE(6,*) ' NUMBER OF SPECTRA=',NDET
      WRITE(6,*) ' Spectrum number=',ISPECT
      WRITE(6,*) ' Number of monitor spectrum=',nmon 
* DETERMINE WHICH SPECTRA IS OUT,THIN,THICK
      IOUT=1
      ITHIN=2
      ITHICK=3

       DO IP=1,NPD
        DO IT=2,NCHAN
         IF(T(IT).GT.100.0.AND.T(IT).LT.400.0) THEN 
          SUM1(IP)=SUM1(IP)+ICTA(IP,IT)
         ELSE IF(T(IT).GT.550.0.AND.T(IT).LT.600.0) THEN
          SUM2(IP)=SUM2(IP)+ICTA(IP,IT)
         END IF
        END DO
        write(6,*) ' sum1=',sum1(ip),' sum2=',sum2(ip)
       END DO 
       
       DO IP=1,NPD
        NI(IP)=IP
        IF(SUM2(IP).NE.0.0) THEN
         SUM1(IP)=SUM1(IP)/SUM2(IP)
         SUM4(IP)=SUM1(IP)
        END IF
       END DO

       CALL ORDER(SUM1,NI,NPD)
       CALL WRITESP(NI,ISPECT,NDET,NPD) ! Write, out,thin,thick spectra.

*      WRITE(6,*) ' 1'
* Calculate total foil out, thin, thick and convert to Cts/micro-sec        

       IF(ISPECT.GE.135.AND.ISPECT.LE.142) THEN
       NI(1)=ISPECT
       NI(2)=ISPECT+2.0*(NDET+1)
       NI(3)=ISPECT+4.0*(NDET+1)
       NI(4)=ISPECT+1.0*(NDET+1)
       NI(5)=ISPECT+3.0*(NDET+1)
       NI(6)=ISPECT+5.0*(NDET+1)
       ELSE IF(ISPECT.GE.143.AND.ISPECT.LE.150) THEN
       NI(1)=ISPECT+1.0*(NDET+1)
       NI(2)=ISPECT+3.0*(NDET+1)
       NI(3)=ISPECT+5.0*(NDET+1)
       NI(4)=ISPECT
       NI(5)=ISPECT+2.0*(NDET+1)
       NI(6)=ISPECT+4.0*(NDET+1)
       ELSE IF(ISPECT.GE.151.AND.ISPECT.LE.158) THEN
       NI(1)=ISPECT
       NI(2)=ISPECT+2.0*(NDET+1)
       NI(3)=ISPECT+4.0*(NDET+1)
       NI(4)=ISPECT+1.0*(NDET+1)
       NI(5)=ISPECT+3.0*(NDET+1)
       NI(6)=ISPECT+5.0*(NDET+1)
       ELSE IF(ISPECT.GE.159.AND.ISPECT.LE.166) THEN
       NI(1)=ISPECT+1.0*(NDET+1)
       NI(2)=ISPECT+3.0*(NDET+1)
       NI(3)=ISPECT+5.0*(NDET+1)
       NI(4)=ISPECT
       NI(5)=ISPECT+2.0*(NDET+1)
       NI(6)=ISPECT+4.0*(NDET+1)
       ELSE IF(ISPECT.GE.167.AND.ISPECT.LE.174) THEN
       NI(1)=ISPECT
       NI(2)=ISPECT+2.0*(NDET+1)
       NI(3)=ISPECT+4.0*(NDET+1)
       NI(4)=ISPECT+1.0*(NDET+1)
       NI(5)=ISPECT+3.0*(NDET+1)
       NI(6)=ISPECT+5.0*(NDET+1)
       ELSE IF(ISPECT.GE.175.AND.ISPECT.LE.182) THEN
       NI(1)=ISPECT+1.0*(NDET+1)
       NI(2)=ISPECT+3.0*(NDET+1)
       NI(3)=ISPECT+5.0*(NDET+1)
       NI(4)=ISPECT
       NI(5)=ISPECT+2.0*(NDET+1)
       NI(6)=ISPECT+4.0*(NDET+1)
       ELSE IF(ISPECT.GE.183.AND.ISPECT.LE.190) THEN
       NI(1)=ISPECT
       NI(2)=ISPECT+2.0*(NDET+1)
       NI(3)=ISPECT+4.0*(NDET+1)
       NI(4)=ISPECT+1.0*(NDET+1)
       NI(5)=ISPECT+3.0*(NDET+1)
       NI(6)=ISPECT+5.0*(NDET+1)
       ELSE IF(ISPECT.GE.191.AND.ISPECT.LE.198) THEN
       NI(1)=ISPECT+1.0*(NDET+1)
       NI(2)=ISPECT+3.0*(NDET+1)
       NI(3)=ISPECT+5.0*(NDET+1)
       NI(4)=ISPECT
       NI(5)=ISPECT+2.0*(NDET+1)
       NI(6)=ISPECT+4.0*(NDET+1)
       END IF
       WRITE(6,*) ' FOIL OUT SPECTRA ARE'
       WRITE(6,*)  NI(4),NI(5),NI(6)
       WRITE(6,*) ' FOIL IN SPECTRA ARE'
       WRITE(6,*) NI(1),NI(2),NI(3)


       IF(ISPECT.GE.135.AND.ISPECT.LE.142) THEN
       NI(1)=1
       NI(2)=3
       NI(3)=5
       NI(4)=2
       NI(5)=4
       NI(6)=6
       ELSE IF(ISPECT.GE.143.AND.ISPECT.LE.150) THEN
       NI(1)=2
       NI(2)=4
       NI(3)=6
       NI(4)=1
       NI(5)=3
       NI(6)=5
       ELSE IF(ISPECT.GE.151.AND.ISPECT.LE.158) THEN
       NI(1)=1
       NI(2)=3
       NI(3)=5
       NI(4)=2
       NI(5)=4
       NI(6)=6
       ELSE IF(ISPECT.GE.159.AND.ISPECT.LE.166) THEN
       NI(1)=2
       NI(2)=4
       NI(3)=6
       NI(4)=1
       NI(5)=3
       NI(6)=5
       ELSE IF(ISPECT.GE.167.AND.ISPECT.LE.174) THEN
       NI(1)=1
       NI(2)=3
       NI(3)=5
       NI(4)=2
       NI(5)=4
       NI(6)=6
       ELSE IF(ISPECT.GE.175.AND.ISPECT.LE.182) THEN
       NI(1)=2
       NI(2)=4
       NI(3)=6
       NI(4)=1
       NI(5)=3
       NI(6)=5
       ELSE IF(ISPECT.GE.183.AND.ISPECT.LE.190) THEN
       NI(1)=1
       NI(2)=3
       NI(3)=5
       NI(4)=2
       NI(5)=4
       NI(6)=6
       ELSE IF(ISPECT.GE.191.AND.ISPECT.LE.198) THEN
       NI(1)=2
       NI(2)=4
       NI(3)=6
       NI(4)=1
       NI(5)=3
       NI(6)=5
       END IF
!       WRITE(6,*) ' FOIL OUT SPECTRA ARE'
!       WRITE(6,*)  NI(4),NI(5),NI(6)
!       WRITE(6,*) ' FOIL IN SPECTRA ARE'
!       WRITE(6,*) NI(1),NI(2),NI(3)

       DO I=2,NCHAN
        COUT(I)=FLOAT(ICTA(NI(4),I)+ICTA(NI(5),I)+ICTA(NI(6),I))
     $/(T(I)-T(I-1))
        ERROUT(I)=SQRT(FLOAT(ICTA(NI(4),I)+ICTA(NI(5),I)+ICTA(NI(6),I)))
     $/(T(I)-T(I-1))
        CTHIN(I)=FLOAT(ICTA(NI(1),I)+ICTA(NI(2),I)+ICTA(NI(3),I)
     $)/(T(I)-T(I-1))
        ERRTHIN(I)=SQRT(FLOAT(ICTA(NI(1),I)+ICTA(NI(2),I)
     $ +ICTA(NI(3),I)))/(T(I)-T(I-1))
        CTHICK(I)=FLOAT(ICTA(NI(1),I)+ICTA(NI(2),I)
     $)/(T(I)-T(I-1))
        ERRTHICK(I)=SQRT(FLOAT(ICTA(NI(1),I)+ICTA(NI(2),I)
     $))/(T(I)-T(I-1))
       END DO
      SUM3(IOUT)=SUM2(NI(4))+SUM2(NI(5))+SUM2(NI(6))
      SUM3(ITHIN)=SUM2(NI(1))+SUM2(NI(2))+SUM2(NI(3))
      SUM3(ITHICK)=SUM2(NI(1))+SUM2(NI(2))

* Calculate total foil out, thin, thick monitor cts         
      IF(NPD.EQ.2) THEN
       DO I=2,NCHAN
        MOUT(I)=FLOAT(MICTA(NI(2),I))/(T(I)-T(I-1))
        MTHIN(I)=FLOAT(MICTA(NI(1),I))/(T(I)-T(I-1))
       END DO
      ELSE IF(NPD.EQ.3) THEN
       DO I=2,NCHAN
        MOUT(I)=FLOAT(MICTA(NI(3),I))/(T(I)-T(I-1))
        MTHIN(I)=FLOAT(MICTA(NI(2),I))/(T(I)-T(I-1))
        MTHICK(I)=FLOAT(MICTA(NI(1),I))/(T(I)-T(I-1))
       END DO
      ELSE IF(NPD.EQ.6) THEN
       DO I=2,NCHAN
        MOUT(I)=FLOAT(MICTA(NI(5),I)+MICTA(NI(6),I)
     $)/(T(I)-T(I-1))
        MTHIN(I)=FLOAT(MICTA(NI(3),I)+MICTA(NI(4),I)
     $)/(T(I)-T(I-1))
        MTHICK(I)=FLOAT(MICTA(NI(1),I)+MICTA(NI(2),I)
     $)/(T(I)-T(I-1))
       END DO
      ELSE IF(NPD.EQ.9) THEN
       DO I=2,NCHAN
        MOUT(I)=FLOAT(MICTA(NI(7),I)+MICTA(NI(8),I)
     $   +MICTA(NI(9),I))/(T(I)-T(I-1))
        MTHIN(I)=FLOAT(MICTA(NI(4),I)+MICTA(NI(5),I)
     $   +MICTA(NI(6),I))/(T(I)-T(I-1))
        MTHICK(I)=FLOAT(MICTA(NI(1),I)+MICTA(NI(2),I)
     $   +MICTA(NI(3),I))/(T(I)-T(I-1))
       END DO
      END IF     


* Calculate number of frames in each period.
      SUM3T=0.0
      DO I=1,3
       SUM3T=SUM3T+SUM3(I)
      END DO
      DO I=1,3
       if(sum3t.ne.0.0)NFR(I)=NGOOD_FRAMES*SUM3(I)/SUM3T
       IF(NFR(I).EQ.0.0) THEN
         write(6,*) ' nframes=0 in period',i
         nfr(i)=1
       end if
      END DO
      

* Start of routines for dead-time correction
* CONVERT TO COUNTS/FRAME
*      DO I=2,NCHAN
*       COUT(I)=COUT(I)/NFR(IOUT)
*       ERROUT(I)=ERROUT(I)/NFR(IOUT)
*       CTHIN(I)=CTHIN(I)/NFR(ITHIN)
*       ERRTHIN(I)=ERRTHIN(I)/NFR(ITHIN)
*       IF(NPD.NE.2) THEN
*        CTHICK(I)=CTHICK(I)/NFR(ITHICK)
*        ERRTHICK(I)=ERRTHICK(I)/NFR(ITHICK)
*       END IF
*      END DO
*
* Correct for dead times in detector
*      TAU=1.0
*      DO I=2,NCHAN
*       COUT(I)=COUT(I)/(1.0-TAU*COUT(I))
*       ERROUT(I)=ERROUT(I)/(1.0-TAU*ERROUT(I))
*       CTHIN(I)=CTHIN(I)/(1.0-TAU*CTHIN(I))
*       ERRTHIN(I)=ERRTHIN(I)/(1.0-TAU*ERRTHIN(I))
*       IF(NPD.NE.2) THEN
*        ERRTHICK(I)=ERRTHICK(I)/(1.0-TAU*ERRTHICK(I))
*       END IF
*      END DO
*
* CONVERT BACK TO TOTAL COUNTS.
*     DO I=2,NCHAN
*       COUT(I)=COUT(I)*NFR(IOUT)
*       ERROUT(I)=ERROUT(I)*NFR(IOUT)
*       CTHIN(I)=CTHIN(I)*NFR(ITHIN)
*       ERRTHIN(I)=ERRTHIN(I)*NFR(ITHIN)
*       IF(NPD.NE.2) THEN
*        CTHICK(I)=CTHICK(I)*NFR(ITHICK)
*        ERRTHICK(I)=ERRTHICK(I)*NFR(ITHICK)
*       END IF
*      END DO



* Calculate monitor cts between 600 and 700 usec.       
      RMOUT=0.0
      RMTHIN=0.0
      RMTHICK=0.0
      DO I=2,NCHAN
       IF(T(I).GE.600.0.AND.T(I).LT.700.0) THEN
        RMOUT=RMOUT+MOUT(I)
        RMTHIN=RMTHIN+MTHIN(I)
        RMTHICK=RMTHICK+MTHICK(I)
*        WRITE(6,*) ' T=',T(I),' MOUT=',MOUT(I)
       END IF
      END DO
      if(rmthick.eq.0.0)rmthick=1e-6

*      WRITE(6,*) ' RMOUT=',RMOUT,' RMTHIN=',RMTHIN,' RMTHICK=',RMTHICK

* Normalise spectra to monitor counter.
      const=1000.0 ! Constant to give suitable scale in tof spectra.
      DO I=2,NCHAN
       COUT(I)=COUT(I)*const/RMOUT
       CTHIN(I)=CTHIN(I)*CONST/RMTHIN
       CTHICK(I)=CTHICK(I)*CONST/RMTHICK
       ERROUT(I)=ERROUT(I)*CONST/RMOUT
       ERRTHIN(I)=ERRTHIN(I)*CONST/RMTHIN
       ERRTHICK(I)=ERRTHICK(I)*CONST/RMTHICK
      END DO


* Calculate cts between 400 and 450 usec.       
      SUMOUT=0.0
      SUMTHIN=0.0
      SUMTHICK=0.0
      DO I=2,NCHAN
       IF(T(I).GE.400.0.AND.T(I).LT.450.0) THEN
        SUMOUT=SUMOUT+COUT(I)
        SUMTHIN=SUMTHIN+CTHIN(I)
        SUMTHICK=SUMTHICK+CTHICK(I)
       END IF
      END DO


      if(sumthin.eq.0.0) then
        write(6,*) ' no cts in thin foil spectra'
        sumthin=1.0
      end if
      if(sumout.eq.0.0) then
        write(6,*) ' no cts in foil out spectra'
        sumout=1.0
      end if
      if(sumthick.eq.0.0) then
        write(6,*) ' no cts in thick foil spectra'
        sumthick=1.0
      end if
 
* Normalise spectra to same area as foil out between 400 and 450 usec.
      DO I=2,NCHAN
       CTHIN(I)=CTHIN(I)*SUMOUT/SUMTHIN
       ERRTHIN(I)=ERRTHIN(I)*SUMOUT/SUMTHIN
       CTHICK(I)=CTHICK(I)*SUMOUT/SUMTHICK
       ERRTHICK(I)=ERRTHICK(I)*SUMOUT/SUMTHICK
      END DO

      IT=0.0
      IF(IOPT.EQ.1)THEN ! THIN DIFFERENCE
       DT=T(3)-T(2)
       DO I=2,NCHAN
        IF(T(I).LT.600.0.and.t(i).gt.1.0) THEN
         IT=IT+1
         XOUT(IT)=T(I)
         YOUT(IT)=-COUT(I)+CTHIN(I)
         EOUT(IT)=SQRT(ERROUT(I)**2+ERRTHIN(I)**2)
        END IF
       END DO
      ELSE IF(IOPT.EQ.2) THEN ! DOUBLE DIFFERENCE
       DO I=2,NCHAN
        IF(T(I).LT.600.0) THEN
         IT=IT+1
         XOUT(IT)=T(I)
         YOUT(IT)=COUT(I)*(1-BETA)-CTHIN(I)+BETA*CTHICK(I)
         EOUT(I)=SQRT((1-BETA)**2*ERROUT(I)**2+ERRTHIN(I)**2
     $           +BETA**2*ERRTHICK(I)**2)
        END IF
       END DO
      ELSE IF(IOPT.EQ.3) THEN ! THICK DIFFERENCE
       DO I=2,NCHAN
        IF(T(I).LT.600.0) THEN
         IT=IT+1
         XOUT(IT)=T(I)
         YOUT(IT)=COUT(I)-CTHICK(I)
         EOUT(IT)=SQRT(ERROUT(I)**2+ERRTHICK(I)**2)
       END IF
       END DO
      END IF 
         


      LPTOUT=IT-1
      write(4,*) ispect,lptout
      do i=1,lptout
       if(t(i).lt.600.0) then
        write(4,*) xout(i),yout(i),eout(i)
       end if
      end do 

 1000 Continue
      write(4,*) irun1,irun2
      close(4)
      write(6,*) ' Data is in TEMP.DAT'
      END

c Sort in order of increasing x values.
      SUBROUTINE ORDER(X,Y,N)
      REAL X(N)
      INTEGER Y(N)

      DO IS=1,N
       IMIN=IS
       XMIN=X(IS)
       DO I=IS,N 
        IF(X(I).LE.XMIN) THEN
        IMIN=I
        XMIN=X(I)
        END IF
       END DO
* Put smallest x value and corresponding y value in X(IS) and Y(IS).
       XTEMP=X(IMIN)
       YTEMP=Y(IMIN)
       X(IMIN)=X(IS)
       Y(IMIN)=Y(IS)
       X(IS)=XTEMP
       Y(IS)=YTEMP
      END DO

      END

      SUBROUTINE WRITESP(NI,ISPECT,NDET,NPD)
      INTEGER NI(9)
       If(npd.eq.9) then
       write(6,*) ' Foil out Spectra are:'
       do i=7,9
         IV=ISPECT+(ni(i)-1)*(NDET+1)
         write(6,*) iv
       end do

       write(6,*) ' Thin Foil Spectra are:'
       do i=4,6
         IV=ISPECT+(ni(i)-1)*(NDET+1)
         write(6,*) iv
       end do

       write(6,*) ' Thick foil Spectra are:'
       do i=1,3
         IV=ISPECT+(ni(i)-1)*(NDET+1)
         write(6,*) iv
       end do
      end if

       If(npd.eq.3) then
       write(6,*) ' Foil out Spectrum is:'
         IV=ISPECT+(ni(3)-1)*(NDET+1)
         write(6,*) iv
       write(6,*) ' Thin foil Spectrum is:'
         IV=ISPECT+(ni(2)-1)*(NDET+1)
         write(6,*) iv
       write(6,*) ' Thick foil Spectrum is:'
         IV=ISPECT+(ni(1)-1)*(NDET+1)
         write(6,*) iv
       end if

       If(npd.eq.2) then
       write(6,*) ' Foil out Spectrum is:'
         IV=ISPECT+(ni(2)-1)*(NDET+1)
         write(6,*) iv
       write(6,*) ' Thin foil Spectrum is:'
         IV=ISPECT+(ni(1)-1)*(NDET+1)
         write(6,*) iv
       end if
       RETURN
       END

*  FUNCTION TO READ NUMBER OF MONITOR COUNTER
      function mon(runid)
      integer nmon*4,nrun*4
      CHARACTER runid*80

      CALL GETPARI(runid,'MDET',NMON,100,LOUT,ERRCODE)
      MON=NMON
      END