      program ana_scorj ! Bernd Berg Feb 20 2012.
C Analysis of scalar field correlation functions.
C Prepares plot data on c0.d and c1.d. 
C Generation of nrpt jackknife file for mass fit.
      include '../../Libs/Fortran/implicit.sta'
      character cd*1,clx*3,clt*3,cb*3 ! cb labels bins.
      include '../../Libs/Fortran/constants.par'
      include 'lat.par'
      include 'mc.par'
C For autocorrelation calculations:
C nbins  Number of bins used for jackknife. 
C        Normally nbins=nrpt. Can be diferent.
      parameter(nbins=nrpt,lmean=.true.)
C Set nth<mth if required to CUT OUT NOISE for massfit.
      parameter(nth=mth)
      include '../../Libs/Fortran_par/common_lat.f'
      include '../../Libs/Scalar_par/common_scalar.f'
      dimension scorj0(nrpt,0:mth),scorj1(nrpt,0:mth)
      dimension scor0(nrpt,0:mth),scor1(nrpt,0:mth),eb0(0:mth),
     &  cor1(0:mth),eb1(0:mth),acpt(nrpt)
      include 'lat.dat'
      ltest=.false.
      ltest=.true.
C
      write(iuo,'(/," ana_scorj: nbins,lmean =",I4,L2)') nbins,lmean
      write(iuo,'(/," Scalar Action")')
      write(cd,'(I1.1)') nd ! Prepare to read the data.
      write(clx,'(I3.3)') nla(1)
      write(clt,'(I3.3)') nla(nd)
      write(cfile,'(A,"d",A,"x",A,"t")') cd,clx,clt
      open(iud1,file=cfln//cfile//".D",form='unformatted',status="old")
      read(iud1) xm2_in,xla_in,nd_in,nla,nlink,nequi_in,nrpt_in,
     &           nmeas_in,nsw_in
      write(iuo,'(/," xm2,xla,nlink,nd:",2F15.10,I10,I3,2X)')
     &                xm2,xla,nlink,nd_in
      write(iuo,'(/," nla:",4I5)') nla
      write(iuo,'(" nequi,nrpt,nmeas,nsw: ",I15,2I10,I6)')
     &              nequi_in,nrpt_in,nmeas_in,nsw_in
      if(nd.ne.nd_in) stop "nd.ne.nd_in."
      if(nequi.ne.nequi_in) stop "nequi.ne.nequi_in."
      if(nrpt.ne.nrpt_in) stop "nrpt.ne.nrpt_in."
      if(nmeas.ne.nmeas_in) stop "nmeas.ne.nmeas_in."
      if(nsw.ne.nsw_in) stop "nsw.ne.nsw_in."
C
      do irpt=1,nrpt
        read(iud1) tsann,tsa2,tsphi,acpt(irpt),phi0sum,
     &             cor0,cor1r,cor1i,i_in
        do idel=0,mth
          scor0(irpt,idel)=cor0(idel)
          scor1(irpt,idel)=half*(cor1r(idel)+cor1i(idel))
        enddo
        if(irpt.ne.i_in) stop "irpt.ne.i_in incomplete read."
      enddo
      close(iud1)
C
c     call steb0(nrpt,act,actm,av,actme)
c     write(iuo,'(" actm = ",F15.10," +/-",F15.10)') actm,actme
      call steb0(nrpt,acpt,acptm,av,acptme)
      write(iuo,'(" acptm =",F15.10," +/-",F15.10)') acptm,acptme
C      
      do idel=0,mth
        call steb0(nrpt,scor0(1,idel),cor0(idel),v0,eb0(idel))
        call steb0(nrpt,scor1(1,idel),cor1(idel),v0,eb1(idel))
        call datjack(nrpt,scor0(1,idel),scorj0(1,idel))
        call datjack(nrpt,scor1(1,idel),scorj1(1,idel))
      enddo
      open(iud1,file="c0j000.d",form='formatted',status='unknown')
      open(iud2,file="c1j000.d",form='formatted',status='unknown')
      open(iud3,file="c0.d",form='formatted',status='unknown')
      open(iud4,file="c1.d",form='formatted',status='unknown')
      write(iud1,'(I10)') 0
      xmass0=log(cor0(0)/cor0(1))
      write(iud1,'(I10,4X,2G14.6)') (nth+1),cor0(0),xmass0
      write(iud2,'(I10)') 1
      xmass1=log(cor1(0)/cor1(1))
      write(iud2,'(I10,4X,2G14.6)') (nth+1),cor1(0),xmass1
      do idel=0,nth ! write only up to nth.
        xdel=idel*one
        write(iud1,'(3G14.6)') xdel,cor0(idel),eb0(idel)
        write(iud2,'(3G14.6)') xdel,cor1(idel),eb1(idel)
        write(iud3,'(3G14.6)') xdel,cor0(idel),eb0(idel)
        write(iud4,'(3G14.6)') xdel,cor1(idel),eb1(idel)
      enddo
      close(iud1)
      close(iud2)
      close(iud3)
      close(iud4)
      if(ltest) stop "ana_scorj: ltest."
C
      do irpt=1,nrpt ! Write jackknife data files.
       write(cb,'(I3.3)') irpt
       open(iud1,file="c0j"//cb//".d",form='formatted',status='unknown')
       open(iud2,file="c1j"//cb//".d",form='formatted',status='unknown')
       write(iud1,'(I10)') 0
       write(iud1,'(I10,4X,2G14.6)') (nth+1),cor0(0),xmass0
       write(iud2,'(I10)') 1
       write(iud2,'(I10,4X,2G14.6)') (nth+1),cor1(0),xmass1
       do idel=0,nth ! write only up to nth.
         xdel=idel*one
         write(iud1,'(3G14.6)') xdel,scorj0(irpt,idel),eb0(idel)
         write(iud2,'(3G14.6)') xdel,scorj1(irpt,idel),eb1(idel)
       enddo
      enddo
C
      stop
      end

      include '../../Libs/Fortran/datjack.f'
      include '../../Libs/Fortran/steb0.f'
      include '../../Libs/Fortran/stmean.f'
