* This program is used as an example for rebinning a histograms of say 200 
* bins to a histogram of 1500 bins
* Using a text file containing the first histogram, interpolating a 
* function out of it and then putting that function into a histograms of 
* 1500 bins.
* NEED TO BE CAREFUL TO THE NORMALIZATION -- Check with locate -i
* I want to go from 75MeV bins to 10 MeV bins
* Need to divide my output by 7.5
*
* It is good practice to always check the overall normalization of the 
* input and output files and make sure they agree.
* Also you might need to change the order of the polynomial that you used
* so plot the text file and the output of DIVDIF before rebinning to check
* that it looks fine.



      program crs_txt_hbk
      IMPLICIT NONE
      integer nwpawc
     
      parameter (nwpawc=100000)
      integer h
      common /PAWC/h(nwpawc)
 

      integer i
      real enu(200),crsnue(200),crsnueb(200),flux(200)
      real EBINCENTER,XSEC,EBINWIDTH,divdif,flux_out
      
      call hlimit(NWPAWC)
 
      call hbook1(201,'crs(nue) QE       ',1500, 0.0, 15.0, 0.0)
      call hbook1(202,'crs(anti-nue) QE  ',1500, 0.0, 15.0, 0.0)
      call hbook1(400,'Flux',1500, 0.0, 15.0, 0.0)
      
      open(90,file='crsqe.txt',form='formatted',status='old')
      open(91,file='28gev_numi.txt',form='formatted',status='old')

      do i=1,200
         enu(i)=0
         crsnue(i)=0
         crsnueb(i)=0
         flux(i)=0
         read(90,*)enu(i),crsnue(i),crsnueb(i)
         read(91,*)flux(i)
      enddo

      close (90)
      close (91)
      ebinwidth=0.010  ! final binwidth that I want

c
c      Use polynomial interpolation to fill histogram using table values:
c      ==================================================================
       DO I = 1, 1500
          xsec=0
          flux_out=0
          EBINCENTER = EBINWIDTH*(FLOAT(I)-0.5)
          write(*,*) i,ebincenter,ebinwidth

c         Flux files:
c         -------------
          flux_out= DIVDIF(flux,enu,200,ebincenter,10)
          CALL HF1(400,EBINCENTER,flux_out/7.5)



c
c         Neutrino cross section
c         ----------------------
          xsec= DIVDIF(crsnue,enu,200,ebincenter,10)
          CALL HF1(201,EBINCENTER,XSEC/7.5)

          write(*,*) ebincenter,xsec/7.5

c

          xsec=0
c         Antineutrino cross section
c         --------------------------
          xsec= DIVDIF(crsnueb,enu,200,ebincenter,10)
          CALL HF1(202,EBINCENTER,XSEC/7.5)
          write(*,*) ebincenter,xsec/7.5

       END DO



       call hrput(0,'crsqe.hbk','NT')

       end