program miniboone_exampleoscfit_nubar c Description: this is an example FORTRAN program to illustrate usage c of MiniBooNE electron and muon antineutrino public data from May 2009 c data release, to perform a 2-neutrino muon-to-electron antineutrino c oscillation fit c Authors: MiniBooNE collaboration c Date: May 2009 implicit none c ------------------------------------------------------------------ c Declarations below c number of enuqe (reconstructed electron (anti)neutrino energy) bins integer nbins_nue integer nbins_numu parameter (nbins_nue=8) parameter (nbins_numu=8) c c 1D array of enuqe bin boundaries for electron antineutrino event c distribution (MeV) double precision enuqevec(nbins_nue+1) c c 1D array of measured nue events per enuqe bin integer nuedata(nbins_nue) c c 1D array of measured numu CCQE events per enuqe bin integer numudata(nbins_numu) c c 1D side-by-side array of measured nue events and measured c numu CCQE events per enuqe bin integer fulldata(nbins_nue+nbins_numu) c c 1D array of predicted nue background per enuqe bin double precision nuebgr(nbins_nue) c c 1D array of predicted numu CCQE events per enuqe bin double precision numu(nbins_numu) c c 1D side-by-side array of predicted numu->nue oscillation and c nue background events, and numu CCQE events per enuqe bin double precision prediction(nbins_nue+nbins_numu) 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 (anti)neutrino energy (MeV) c enutrue: true (anti)neutrino energy (MeV) c lnutrue: distance between (anti)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=28685) 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_nue) 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_nue) c 2D array of 3x3-block fractional covariance matrix for: full numu->nue c oscillation enuqe distribution (systematic errors only), c nue background enuqe distribution (systematic + statistical c errors), and numu CCQE enuqe distribution (systematic + c statistical errors) c (This corresponds to the three distributions, side-by-side) double precision full_fractcovmatrix(nbins_nue+nbins_nue+ + nbins_numu,nbins_nue+nbins_nue+nbins_numu) c 2D array of corresponding full covariance matrix, scaled to c signal prediction double precision full_covmatrix(nbins_nue+nbins_nue+nbins_numu, + nbins_nue+nbins_nue+nbins_numu) c c covariance matrix and its inverse, used in chi2 calculation, c rebinned double precision covmatrix(nbins_nue+nbins_numu, + nbins_nue+nbins_numu) double precision inversecovmatrix(nbins_nue+nbins_numu, + nbins_nue+nbins_numu) 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/nuesignal_unitsinsq2theta, + prediction,inversecovmatrix, + chi2null,chi2boundary,sensitivity, + fulldata 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=190) double precision log10dm2min,log10dm2max parameter (log10dm2min=-2.0d0,log10dm2max=2.0d0) 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_nubar.txt', + status='old') read(10,*) (enuqevec(i),i=1,nbins_nue+1) c c get measured nue events per enuqe bin open(11,file='data/default/miniboone_nuebardata.txt', + status='old') read(11,*) (nuedata(i),i=1,nbins_nue) c c get measured numu CCQE events per enuqe bin open(12,file='data/default/miniboone_numubardata.txt', + status='old') read(12,*) (numudata(i),i=1,nbins_numu) c c get predicted nue background events per enuqe bin open(13,file='data/default/miniboone_nuebarbgr.txt', + status='old') read(13,*) (nuebgr(i),i=1,nbins_nue) c c get predicted numu CCQE events per enuqe bin open(14,file='data/default/miniboone_numubar.txt', + status='old') read(14,*) (numu(i),i=1,nbins_numu) c c get fractional covariance matrix for full (100%) numu->nue c oscillation events (sys-only), and nue background and numu c CCQE enuqe distributions (stat + syst) open(15,file= + 'data/default/miniboone_full_fractcovmatrix_nubar.txt', + status='old') do i=1,nbins_nue+nbins_nue+nbins_numu read(15,*) (full_fractcovmatrix(i,j), + j=1,nbins_nue+nbins_nue+nbins_numu) enddo c c get nfulloscevts full (100%) numu->nue oscillation events c after nue cuts open(16,file= + 'data/default/miniboone_numubarnuebarfullosc_ntuple.txt', + status='old') do ifulloscevt=1,nfulloscevts read(16,*) numunuefullosc_enuqe(ifulloscevt), + numunuefullosc_enutrue(ifulloscevt), + numunuefullosc_lnutrue(ifulloscevt), + numunuefullosc_weight(ifulloscevt) enddo c c open file where to store dm2,sinsq2thetamax values open(17,file= + 'data/default/miniboone_limit_sensitivity_90cl_nubar.txt', + status='unknown') chi2null=1.0d7 chi2boundary=1.0d7 c do i=1,nbins_nue fulldata(i)=nuedata(i) prediction(i)=nuebgr(i) enddo do i=1,nbins_numu fulldata(i+nbins_nue)=numudata(i) prediction(i+nbins_nue)=numu(i) enddo do i=1,nbins_nue+nbins_nue+nbins_numu do j=1,nbins_nue+nbins_nue+nbins_numu full_covmatrix(i,j)=0.d0 enddo enddo 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_nue nuesignal_unitsinsq2theta(i)=0.d0 enddo c do ifulloscevt=1,nfulloscevts do i=1,nbins_nue 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_nue 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,3.0d-4,1.0d0,ierflag) do i=1, nbins_nue+nbins_nue+nbins_numu do j=1, nbins_nue+nbins_nue+nbins_numu full_covmatrix(i,j)=0.d0 enddo enddo do i=nbins_nue, nbins_nue+nbins_nue+nbins_numu do j=nbins_nue, nbins_nue+nbins_nue+nbins_numu full_covmatrix(i,j)=full_fractcovmatrix(i,j)* + prediction(i-nbins_nue)* + prediction(j-nbins_nue) enddo enddo do i=1, nbins_nue+nbins_numu do j=1, nbins_nue+nbins_numu if (i.le.nbins_nue .and. j.le.nbins_nue) then covmatrix(i,j)=full_covmatrix(i,j)+ + full_covmatrix(i+nbins_nue,j)+ + full_covmatrix(i,nbins_nue+j)+ + full_covmatrix(i+nbins_nue,j+nbins_nue) elseif (i.le.nbins_nue .and. j.gt.nbins_nue) then covmatrix(i,j)=full_covmatrix(i,j+nbins_nue)+ + full_covmatrix(i+nbins_nue,j+nbins_nue) elseif (i.gt.nbins_nue .and. j.le.nbins_nue) then covmatrix(i,j)=full_covmatrix(i+nbins_nue,j)+ + full_covmatrix(i+nbins_nue,j+nbins_nue) elseif (i.gt.nbins_nue .and. j.gt.nbins_nue) then covmatrix(i,j)=full_covmatrix(i+nbins_nue,j+nbins_nue) 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 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,3.0d-4,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_nue nuesignal_bestfitsinsq2theta(i)=0.d0 enddo chi2min(0)=1.d7 hasconverged=.false. do iiter=1,nitersmax c build covariance matrix do i=1,nbins_nue+nbins_nue+nbins_numu do j=1,nbins_nue+nbins_nue+nbins_numu if (i.le.nbins_nue .and. j.le.nbins_nue) then full_covmatrix(i,j)=full_fractcovmatrix(i,j)* + nuesignal_bestfitsinsq2theta(i)* + nuesignal_bestfitsinsq2theta(j) elseif (i.le.nbins_nue .and. j.gt.nbins_nue) then full_covmatrix(i,j)=full_fractcovmatrix(i,j)* + nuesignal_bestfitsinsq2theta(i)* + prediction(j-nbins_nue) elseif (i.gt.nbins_nue .and. j.le.nbins_nue) then full_covmatrix(i,j)=full_fractcovmatrix(i,j)* + prediction(i-nbins_nue)* + nuesignal_bestfitsinsq2theta(j) elseif (i.gt.nbins_nue .and. j.gt.nbins_nue) then full_covmatrix(i,j)=full_fractcovmatrix(i,j)* + prediction(i-nbins_nue)* + prediction(j-nbins_nue) endif c Note: add statistical error of signal prediction if (i.le.nbins_nue .and. j.le.nbins_nue) then if (i.eq.j) then full_covmatrix(i,j)=full_covmatrix(i,j)+ + nuesignal_bestfitsinsq2theta(i) endif endif enddo enddo c do i=1,nbins_nue+nbins_numu do j=1,nbins_nue+nbins_numu if (i.le.nbins_nue .and. j.le.nbins_nue) then covmatrix(i,j)=full_covmatrix(i,j)+ + full_covmatrix(i+nbins_nue,j)+ + full_covmatrix(i,nbins_nue+j)+ + full_covmatrix(i+nbins_nue,j+nbins_nue) elseif (i.le.nbins_nue .and. j.gt.nbins_nue) then covmatrix(i,j)=full_covmatrix(i,j+nbins_nue)+ + full_covmatrix(i+nbins_nue,j+nbins_nue) elseif (i.gt.nbins_nue .and. j.le.nbins_nue) then covmatrix(i,j)=full_covmatrix(i+nbins_nue,j)+ + full_covmatrix(i+nbins_nue,j+nbins_nue) elseif (i.gt.nbins_nue .and. j.gt.nbins_nue) then covmatrix(i,j)=full_covmatrix(i+nbins_nue, + j+nbins_nue) 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 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 once, 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_nue 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(17,*) 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) close(17) return end c ================================================================== double precision function oscprob(enutrue,lnutrue,dm2) c Description: compute oscillation probability for a given dm2, c true (anti)neutrino energy and true (anti)neutrino baseline, c 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_nue integer nbins_numu parameter (nbins_nue=8) parameter (nbins_numu=8) integer fulldata(nbins_nue+nbins_numu) double precision prediction(nbins_nue+nbins_numu) double precision nuesignal_unitsinsq2theta(nbins_nue) double precision inversecovmatrix(nbins_nue+nbins_numu, + nbins_nue+nbins_numu) double precision chi2null,chi2boundary logical sensitivity common /miniboonedata/nuesignal_unitsinsq2theta, + prediction,inversecovmatrix, + chi2null,chi2boundary,sensitivity, + fulldata c 1D array of predicted nue signal events per enuqe bin double precision nuesignal(nbins_nue) 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_nue nuesignal(i)=sinsq2theta*nuesignal_unitsinsq2theta(i) enddo c now form chi2 chi2=0.d0 if (sensitivity) then do i=1,nbins_nue+nbins_numu do j=1,nbins_nue+nbins_numu if (i.le.nbins_nue .and. j.le.nbins_nue) then chi2=chi2+ + (prediction(i)- + (prediction(i)+nuesignal(i)))* + inversecovmatrix(i,j)* + (prediction(j)- + (prediction(j)+nuesignal(j))) elseif (i.le.nbins_nue .and. j.gt.nbins_nue) then chi2=chi2+ + (prediction(i)- + (prediction(i)+nuesignal(i)))* + inversecovmatrix(i,j)* + (prediction(j)- + (prediction(j))) elseif (i.gt.nbins_nue .and. j.le.nbins_nue) then chi2=chi2+ + (prediction(i)- + (prediction(i)))* + inversecovmatrix(i,j)* + (prediction(j)- + (prediction(j)+nuesignal(j))) elseif (i.gt.nbins_nue .and. j.gt.nbins_nue) then chi2=chi2+ + (prediction(i)- + (prediction(i)))* + inversecovmatrix(i,j)* + (prediction(j)- + (prediction(j))) endif enddo enddo else do i=1,nbins_nue+nbins_numu do j=1,nbins_nue+nbins_numu if (i.le.nbins_nue .and. j.le.nbins_nue) then chi2=chi2+ + (dble(fulldata(i))- + (prediction(i)+nuesignal(i)))* + inversecovmatrix(i,j)* + (dble(fulldata(j))- + (prediction(j)+nuesignal(j))) elseif (i.le.nbins_nue .and. j.gt.nbins_nue) then chi2=chi2+ + (dble(fulldata(i))- + (prediction(i)+nuesignal(i)))* + inversecovmatrix(i,j)* + (dble(fulldata(j))- + (prediction(j))) elseif (i.gt.nbins_nue .and. j.le.nbins_nue) then chi2=chi2+ + (dble(fulldata(i))- + (prediction(i)))* + inversecovmatrix(i,j)* + (dble(fulldata(j))- + (prediction(j)+nuesignal(j))) elseif (i.gt.nbins_nue .and. j.gt.nbins_nue) then chi2=chi2+ + (dble(fulldata(i))- + (prediction(i)))* + inversecovmatrix(i,j)* + (dble(fulldata(j))- + (prediction(j))) endif 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_nue integer nbins_numu parameter (nbins_nue=8) parameter (nbins_numu=8) integer i,j double precision matrix(nbins_nue+nbins_numu,nbins_nue+nbins_numu) double precision matrix_tmp(nbins_nue+nbins_numu, + nbins_nue+nbins_numu) c LAPACK variables double precision wr(nbins_nue+nbins_numu) double precision work(3*(nbins_nue+nbins_numu)) integer ok c ------------------------------------------------------------------ c Executable statements below do i=1,nbins_nue+nbins_numu do j=1,nbins_nue+nbins_numu matrix_tmp(i,j)=matrix(i,j) enddo enddo call DSYEV('N','U',nbins_nue+nbins_numu,matrix_tmp, + nbins_nue+nbins_numu,wr,work, + 3*(nbins_nue+nbins_numu),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_nue integer nbins_numu parameter (nbins_nue=8) parameter (nbins_numu=8) integer i,j double precision matrix(nbins_nue+nbins_numu,nbins_nue+nbins_numu) double precision inversematrix(nbins_nue+nbins_numu, + nbins_nue+nbins_numu) c LAPACK variables integer ifail c ------------------------------------------------------------------ c Executable statements below do i=1,nbins_nue+nbins_numu do j=1,nbins_nue+nbins_numu inversematrix(i,j)=matrix(i,j) enddo enddo call DSINV(nbins_nue+nbins_numu,inversematrix, + nbins_nue+nbins_numu,ifail) if (ifail.ne.0) then print *,'Matrix Inversion problems!', + ' dsinv ifail: ',ifail stop endif return end