program miniboone_exampleoscfit
c     Description: this is an example FORTRAN program to illustrate usage
c     of MiniBooNE electron neutrino public information from 2007-04-11
c     data release, to perform a 2-neutrino muon-to-electron neutrino 
c     oscillation fit
c     Authors: MiniBooNE collaboration
c     Date: 2007-04-11
      implicit none

c     ------------------------------------------------------------------

c     Declarations below

c     number of enuqe (reconstructed electron neutrino energy) bins
      integer nbins
      parameter (nbins=8)
c
c     1D array of enuqe bin boundaries (MeV)
      double precision enuqevec(nbins+1)
c
c     1D array of measured nue events per enuqe bin
      integer nuedata(nbins)
c
c     1D array of predicted nue background per enuqe bin
      double precision nuebgr(nbins)
c
c     2D array of fractional covariance matrix for nue background
c     enuqe distribution (stat + syst)
      double precision nuebgr_fractcovmatrix(nbins,nbins)
c
c     ntuple file of nfulloscevts full (100%) numu->nue oscillation 
c     events after nue cuts, with one row per event, and four columns:
c     enuqe, enutrue, lnutrue, weight
c     enuqe: reconstructed neutrino energy (MeV)
c     enutrue: true neutrino energy (MeV)
c     lnutrue: distance between neutrino production and detection points (cm)
c     weight: event weight, with weight normalization such that the sum of
c     all weights divided by the number of events in the file equals the
c     number of expected full oscillation events across all enuqe
      integer nfulloscevts
      parameter (nfulloscevts=9934)
      double precision numunuefullosc_enuqe(nfulloscevts)
      double precision numunuefullosc_enutrue(nfulloscevts)
      double precision numunuefullosc_lnutrue(nfulloscevts)
      double precision numunuefullosc_weight(nfulloscevts)
c

c     1D array of predicted nue signal events per enuqe bin,
c     for sinsq2theta=1
      double precision nuesignal_unitsinsq2theta(nbins)

c     1D array of predicted nue signal events per enuqe bin,
c     at best-fit sinsq2theta for a specific dm2
      double precision nuesignal_bestfitsinsq2theta(nbins)

c     2D array of fractional covariance matrix for full numu->nue
c     oscillation enuqe distribution (systematic errors only)
      double precision numunuefullosc_fractcovmatrix(nbins,nbins)
c
c     covariance matrix and its inverse, used in chi2 calculation
      double precision covmatrix(nbins,nbins)
      double precision inversecovmatrix(nbins,nbins)

c     flag to compute either sensitivity or limit. In the latter case,
c     the sensitivity flag is set to false
      logical sensitivity

c     chi2 for no-oscillations
      double precision chi2null
c     chi2 at sinsq2theta=1 boundary
      double precision chi2boundary

      common /miniboonedata/nuedata,nuebgr,
     +     nuesignal_unitsinsq2theta,inversecovmatrix,
     +     chi2null,chi2boundary,sensitivity

c     dm2 value, and maximum value of sinsq2theta at given confidence
c     level, sinsq2thetamax, and sinsq2theta sensitivity
      double precision dm2,sinsq2thetamax, sinsq2thetasens
c
c     indeces for enuqe bins
      integer i,j
c     index for full numu->nue oscillation events
      integer ifulloscevt
c     
c     oscillation probability function, for sinsq2theta=1
      double precision oscprob

c     logical function to check if covariance matrix and its inverse
c     are positive-definite
      logical matrixispositivedefinite

c     MINUIT variables, for minimization
      double precision arglis(2),futil
      integer ierflag
      external miniboonechi2

c     MINUIT output variables
      double precision fmin,fedm,errdef
      integer npari,nparx,istat
      character*10 chnam
      double precision val,error,bnd1,bnd2
      integer ivarbl
      double precision eplus,eminus,eparab,globcc

c     confidence level and corresponding quantile of chi2 distribution
      real confidenceLevel
      real gausin
      double precision dchi2cut

c     dm2 grid: number of points, dm2_min, dm2_max
      integer ndm2points
      parameter (ndm2points=1000)
      double precision log10dm2min,log10dm2max
      parameter (log10dm2min=-2.0d0,log10dm2max=2.0d0)
c      parameter (log10dm2min=-1.6d0,log10dm2max=-1.2d0)
c     index for dm2 point
      integer idm2point

c     covariance matrix is computed iteratively, using best-fit
c     sinsq2theta from the previous iteration at each dm2 grid point
c     Iteration process ends either after convergence criterion in
c     best-fit chi2 has been reached for all dm2 grid points, as set 
c     by chi2difftolerance, or after niters iterations, whichever is 
c     reached first
c
c     max number of iterations for computing the covariance matrix
      integer nitersmax,iiter
      parameter (nitersmax=5)
c     convergence criterion in best-fit chi2
      double precision chi2difftolerance
      parameter (chi2difftolerance=1.d-2)
      double precision chi2min(0:nitersmax)
      logical hasconverged
      

c     ------------------------------------------------------------------

c     Executable statements below

c     get enuqe bin boundaries information
      open(10,file='data/default/miniboone_binboundaries.txt',
     +     status='old')
      read(10,*) (enuqevec(i),i=1,nbins+1)
c
c     get measured nue events per enuqe bin
      open(11,file='data/default/miniboone_nuedata.txt',
     +     status='old')
      read(11,*) (nuedata(i),i=1,nbins)
c
c     get predicted nue background events per enuqe bin
      open(12,file='data/default/miniboone_nuebgr.txt',
     +     status='old')
      read(12,*) (nuebgr(i),i=1,nbins)
c
c     get fractional covariance matrix for nue background
c     enuqe distribution (stat + syst)
      open(13,file=
     +     'data/default/miniboone_nuebgr_fractcovmatrix.txt',
     +     status='old')
      do i=1,nbins
         read(13,*) (nuebgr_fractcovmatrix(i,j),j=1,nbins)
      enddo
c
c     get nfulloscevts full (100%) numu->nue oscillation events 
c     after nue cuts
      open(14,file=
     +     'data/default/miniboone_numunuefullosc_ntuple.txt',
     +     status='old')
      do ifulloscevt=1,nfulloscevts
         read(14,*) numunuefullosc_enuqe(ifulloscevt),
     +        numunuefullosc_enutrue(ifulloscevt),
     +        numunuefullosc_lnutrue(ifulloscevt),
     +        numunuefullosc_weight(ifulloscevt)
      enddo
c
c     get fractional covariance matrix for full numu->nue
c     oscillation enuqe distribution (stat + syst)
      open(15,file=
     +'data/default/miniboone_numunuefullosc_fractcovmatrix.txt',
     +     status='old')
      do i=1,nbins
         read(15,*) (numunuefullosc_fractcovmatrix(i,j),j=1,nbins)
      enddo

c     open file where to store dm2,sinsq2thetamax values
      open(16,file=
     +     'data/default/miniboone_limit_sensitivity_90cl.txt',
     +     status='unknown')

      chi2null=1.0d7
      chi2boundary=1.0d7

c     Initialize MINUIT
      print *,'Initializing MINUIT...'
      call mninit(5,6,7)
      call mnseti('MiniBooNE oscillation fit')

c     MINUIT settings
      arglis(1) = 2.0
      call mnexcm(miniboonechi2,'SET STRATEGY',arglis,1,ierflag,futil)
c     set minuit error (dchi2 value), for given confidence level and
c     one-sided, one parameter
      confidenceLevel=0.90
      dchi2cut = dble(gausin(confidenceLevel))**2
      arglis(1) = dchi2cut
      call mnexcm(miniboonechi2,'SET ERR',arglis,1,ierflag,futil)


c     start dm2 loop
      do idm2point=1,ndm2points
         dm2=10.d0**(log10dm2min+
     +        (log10dm2max-log10dm2min)*(dble(idm2point)-1.d0)/
     +        (dble(ndm2points)-1.d0))

c     compute number of predicted oscillation signal events per enuqe
c     bin for this dm2 value, assuming sinsq2theta=1.0
c
         do i=1,nbins
            nuesignal_unitsinsq2theta(i)=0.d0
         enddo
c
         do ifulloscevt=1,nfulloscevts
            do i=1,nbins
               if(numunuefullosc_enuqe(ifulloscevt).ge.enuqevec(i).and.
     +              numunuefullosc_enuqe(ifulloscevt).lt.enuqevec(i+1))
     +              then
                  nuesignal_unitsinsq2theta(i)=
     +                 nuesignal_unitsinsq2theta(i)+
     +                 numunuefullosc_weight(ifulloscevt)*
     +                 oscprob(numunuefullosc_enutrue(ifulloscevt),
     +                 numunuefullosc_lnutrue(ifulloscevt),dm2)
                  goto 100
               endif
            enddo
 100        continue
         enddo
c
         do i=1,nbins
            nuesignal_unitsinsq2theta(i)=
     +           nuesignal_unitsinsq2theta(i)/dble(nfulloscevts)
         enddo

c     find sensitivity

c     define MINUIT fit parameters
         call mnexcm(miniboonechi2,'CLEAR',arglis,0,ierflag,futil)
         call mnparm(1,'sinsq2th',0.0d0,1.0d-3,-1.0d0,1.0d0,ierflag)

         do i=1,nbins
            do j=1,nbins
               covmatrix(i,j)=nuebgr_fractcovmatrix(i,j)*nuebgr(i)*
     +              nuebgr(j)
            enddo
         enddo

c     make sure covariance matrix is positive-definite
         if (.not.matrixispositivedefinite(covmatrix)) then
            print *,'covariance matrix has negative eigenvalues!',
     +           ' Exiting'
            stop
         endif
         
c     covariance matrix inversion
         call getinversematrix(covmatrix,inversecovmatrix)
         
c     make sure that also inverse covariance matrix is positive-definite
         if (.not.matrixispositivedefinite(inversecovmatrix)) then
            print *,'inverse covariance matrix has negative',
     +           ' eigenvalues! Exiting'
            stop
         endif

         sensitivity=.true.

c     perform the MINUIT minimization for this dm2 value
         arglis(1) = 10000.d0
         arglis(2) = 1.d00
         call mnexcm(miniboonechi2,'MINIMIZE',arglis,2,ierflag,futil)
         call mnpout(1,chnam,val,error,bnd1,bnd2,ivarbl)

c     Get minimization status
         call mnstat(fmin,fedm,errdef,npari,nparx,istat)

c     Get fit parameter best-fit value
         call mnpout(1,chnam,val,error,bnd1,bnd2,ivarbl)

c     MINUIT MINOS error analysis
         arglis(1) = 10000.d0
         arglis(2) = 1.d0
         call mnexcm(miniboonechi2,'MINOS',arglis,2,ierflag,futil)
 
c     Get fit parameter errors
         call mnerrs(1,eplus,eminus,eparab,globcc)
        
c     Add extra-check to see if sinsq2theta=1 boundary is allowed, 
c     where error estimate with MINOS is difficult
         arglis(1)=1.d0
         arglis(1)=1.d0
         call mnexcm(miniboonechi2,'SET PAR',arglis,2,ierflag,futil)
         call mnexcm(miniboonechi2,'FIX',arglis,1,ierflag,futil) 
         arglis(1)=7.0d0
         call mnexcm(miniboonechi2,'CALL',arglis,1,ierflag,futil)       
         
         sinsq2thetasens=val+eplus
         if (chi2boundary-fmin.lt.dchi2cut) sinsq2thetasens=1.0d0

c     Now fit to actual data, using raster-scan method
         sensitivity=.false.

c     define MINUIT fit parameters
         call mnexcm(miniboonechi2,'CLEAR',arglis,0,ierflag,futil)
         call mnparm(1,'sinsq2th',0.0d0,1.0d-3,-1.0d0,1.0d0,ierflag)

c     begin iterative process, using best-fit sinsq2theta from previous
c     iteration to compute covariance matrix
c     start with no signal events
         do i=1,nbins
            nuesignal_bestfitsinsq2theta(i)=0.d0
         enddo

         chi2min(0)=1.d7
         hasconverged=.false.

         do iiter=1,nitersmax

c     build covariance matrix
            do i=1,nbins
               do j=1,nbins
                  covmatrix(i,j)=nuebgr_fractcovmatrix(i,j)*nuebgr(i)*
     +                 nuebgr(j)+
     +                 numunuefullosc_fractcovmatrix(i,j)*
     +                 nuesignal_bestfitsinsq2theta(i)*
     +                 nuesignal_bestfitsinsq2theta(j)
c     Note: add statistical error of signal prediction
                  if (i.eq.j) then
                     covmatrix(i,j)=covmatrix(i,j)+
     +                    nuesignal_bestfitsinsq2theta(i)
                  endif
               enddo
            enddo

c     make sure covariance matrix is positive-definite
            if (.not.matrixispositivedefinite(covmatrix)) then
               print *,'covariance matrix has negative eigenvalues!',
     +              ' Exiting'
               stop
            endif

c     Covariance matrix inversion.
            call getinversematrix(covmatrix,inversecovmatrix)

c     make sure that also inverse covariance matrix is positive-definite
            if (.not.matrixispositivedefinite(inversecovmatrix)) then
               print *,'inverse covariance matrix has negative',
     +              ' eigenvalues! Exiting'
               stop
            endif

c     compute chi2null, in the case of no oscillation signal events. 
c     The covariance matrix for signal events is the one from the first 
c     iteration
c     Also, need to do this only one, say for first dm2 grid point,
c     since result is dm2-independent
            if (iiter.eq.1.and.idm2point.eq.1) then
               arglis(1)=6.0d0
               call mnexcm(miniboonechi2,'CALL',arglis,1,ierflag,futil)
            endif

c     perform the MINUIT minimization for this dm2 value and covariance
c     matrix
            arglis(1) = 10000.d0
            arglis(2) = 1.d00
            call mnexcm(miniboonechi2,'MINIMIZE',arglis,2,ierflag,futil)
            call mnpout(1,chnam,val,error,bnd1,bnd2,ivarbl)
            
c     compute number of nue signal signal events for best-fit sinsq2theta
            do i=1,nbins
               nuesignal_bestfitsinsq2theta(i)=
     +              val*nuesignal_unitsinsq2theta(i)
            enddo

c     Get minimization status
            call mnstat(fmin,fedm,errdef,npari,nparx,istat)
            chi2min(iiter)=fmin

            if (abs(chi2min(iiter)-chi2min(iiter-1)).lt.
     +           chi2difftolerance) then
               hasconverged=.true.
               goto 200
            endif
                     
         enddo                  ! end of cov matrix calculation iterations

 200     continue

c     Get fit parameter best-fit value
         call mnpout(1,chnam,val,error,bnd1,bnd2,ivarbl)

c     MINUIT MINOS error analysis
         arglis(1) = 100000.d0
         arglis(2) = 1.d0
         call mnexcm(miniboonechi2,'MINOS',arglis,2,ierflag,futil)

c     Get fit parameter errors
         call mnerrs(1,eplus,eminus,eparab,globcc)

c     Add extra-check to see if sinsq2theta=1 boundary is allowed, 
c     where error estimate with MINOS is difficult
         arglis(1)=1.d0
         arglis(1)=1.d0
         call mnexcm(miniboonechi2,'SET PAR',arglis,2,ierflag,futil)
         call mnexcm(miniboonechi2,'FIX',arglis,1,ierflag,futil) 
         arglis(1)=7.0d0
         call mnexcm(miniboonechi2,'CALL',arglis,1,ierflag,futil)       
         
         sinsq2thetamax=val+eplus
         if (chi2boundary-fmin.lt.dchi2cut) sinsq2thetamax=1.0d0

         write(16,*) dm2,sinsq2thetamax,sinsq2thetasens,val,fmin
         
      enddo                     ! end of dm2 loop

      print *, 'chi2null = ',chi2null
      
      close(10)
      close(11)
      close(12)
      close(13)
      close(14)
      close(15)
      close(16)

      return
      end

c     ==================================================================

      double precision function oscprob(enutrue,lnutrue,dm2)
c     Description: compute oscillation probability for a given dm2,
c     true neutrino energy and true neutrino baseline, for sinsq2theta=1
      implicit none

c     ------------------------------------------------------------------

c     Declarations below 

      double precision enutrue,lnutrue,dm2

c     ------------------------------------------------------------------

c     Executable statements below

      if (enutrue.gt.0.) then
c     Note: convert baseline(cm) in baseline(m)
         oscprob=sin(1.267d0*dm2*lnutrue*1.d-2/enutrue)**2
      else
         oscprob=0.5d0
      endif
      return
      end
      

c     ==================================================================

      subroutine miniboonechi2(npar,grad,fval,xval,iflag,futil)
c     Description: miniboonechi2 defines the chi2 function fval being 
c     minimized by minuit
      implicit none

c     ------------------------------------------------------------------

c     Declarations below
      
c     MINUIT variables
      integer npar
c     xval is the vector of (constant and variable) parameters
c     grad is the (optional) vector of first derivatives
      double precision grad(*),xval(*)
c     fval is the calculated function value
c     name of optional utilitary routine
      double precision fval,futil
c     iflag indicated what is to be calculated
      integer iflag
      
c     MiniBooNE data
      integer nbins
      parameter (nbins=8)
      integer nuedata(nbins)
      double precision nuebgr(nbins)
      double precision nuesignal_unitsinsq2theta(nbins)
      double precision inversecovmatrix(nbins,nbins)
      double precision chi2null,chi2boundary
      logical sensitivity

      common /miniboonedata/nuedata,nuebgr,
     +     nuesignal_unitsinsq2theta,inversecovmatrix,
     +     chi2null,chi2boundary,sensitivity

c     1D array of predicted nue signal events per enuqe bin
      double precision nuesignal(nbins)

c     fit parameters
      double precision sinsq2theta

c     indeces for enuqe bins
      integer i,j

c     chi2 function
      double precision chi2

c     ------------------------------------------------------------------

c     Executable statements below

      sinsq2theta=xval(1)

c     number of predicted nue signal events
      do i=1,nbins
         nuesignal(i)=sinsq2theta*nuesignal_unitsinsq2theta(i)
      enddo

c     now form chi2
      chi2=0.d0

      if (sensitivity) then
         do i=1,nbins
            do j=1,nbins
               chi2=chi2+
     +              (dble(nuebgr(i))-(nuebgr(i)+nuesignal(i)))*
     +              inversecovmatrix(i,j)*
     +              (dble(nuebgr(j))-(nuebgr(j)+nuesignal(j)))
            enddo
         enddo
      else        
         do i=1,nbins
            do j=1,nbins
               chi2=chi2+
     +              (dble(nuedata(i))-(nuebgr(i)+nuesignal(i)))*
     +              inversecovmatrix(i,j)*
     +              (dble(nuedata(j))-(nuebgr(j)+nuesignal(j)))
            enddo
         enddo
      endif

      if (iflag.eq.6) chi2null=chi2
      if (iflag.eq.7) chi2boundary=chi2
      fval = chi2

      return
      end

c     ==================================================================

      logical function matrixispositivedefinite(matrix)
c     Description: compute matrix eigenvalues, and check that
c     all eigenvalues are positive, as it should be for a
c     positive-definite matrix
      implicit none

c     ------------------------------------------------------------------

c     Declarations below

      integer nbins
      parameter (nbins=8)
      integer i,j
      double precision matrix(nbins,nbins)
      double precision matrix_tmp(nbins,nbins)

c     LAPACK variables
      double precision wr(nbins)
      double precision work(3*nbins)
      integer ok

c     ------------------------------------------------------------------

c     Executable statements below

      do i=1,nbins
         do j=1,nbins
            matrix_tmp(i,j)=matrix(i,j)
         enddo
      enddo
      call DSYEV('N','U',nbins,matrix_tmp,nbins,wr,work,
     +     3*nbins,ok)
      if (wr(1).lt.0.d0) then
         matrixispositivedefinite=.false.
      else
         matrixispositivedefinite=.true.
      endif

      return
      end

c     ==================================================================

      subroutine getinversematrix(matrix,inversematrix)
c     Description: compute matrix inversion
      implicit none

c     ------------------------------------------------------------------

c     Declarations below

      integer nbins
      parameter (nbins=8)
      integer i,j
      double precision matrix(nbins,nbins)
      double precision inversematrix(nbins,nbins)

c     LAPACK variables
      integer ifail

c     ------------------------------------------------------------------

c     Executable statements below

      do i=1,nbins
         do j=1,nbins
            inversematrix(i,j)=matrix(i,j)
         enddo
      enddo
      call DSINV(nbins,inversematrix,nbins,ifail)
      if (ifail.ne.0) then
         print *,'Matrix Inversion problems!',
     +        ' dsinv ifail: ',ifail
         stop
      endif

      return
      end