subroutine match_covariance_distribution(nmuons,muon_file, muons)

USE nr
use parameters_bmad
use muon_mod
use muon_interface
use random_mod
use sim_utils


IMPLICIT NONE

integer unit, nmuons
character*100 string/'correlation'/
character*120 muon_file
type (muon_struct), allocatable :: muons(:), muons_ref(:)
type (params), allocatable :: gauss(:)
real(rp), allocatable :: p(:), Sigma(:,:), temp(:,:), ave(:)
real(rp), allocatable ::   LT(:,:),  M(:,:)
real(rp), allocatable, save ::  L(:,:)

integer lun, ndim
integer i,j
integer tot/0/, turn/2/
integer nargs
integer nmuons_ref/-1/

logical first/.true./

if(first)then
   
   ndim=8

   allocate(p(1:ndim))
   allocate(Sigma(1:ndim, 1:ndim))
   allocate(L(1:ndim, 1:ndim),LT(1:ndim,1:ndim))
   allocate(M(1:ndim,1:ndim))
   allocate(temp(1:ndim,1:ndim))
   allocate(gauss(1:size(muons)))
   allocate(ave(1:ndim))
   
   lun = lunget()
   open(unit = lun, file = trim(muon_file))

   if(nmuons_ref<0)then
      call get_number_events(lun,  tot)
      nmuons_ref = tot-4
      rewind(unit=lun)
      allocate(muons_ref(1:nmuons_ref))
   endif

   if(.not.allocated(muons_ref))allocate(muons_ref(1:size(muons)))
   call read_phase_space(lun, nmuons_ref, muons_ref, string, tot)

   call compute_correlation(nmuons_ref, muons_ref, Sigma, ave)

   print *,' averge and Covariance matrix'
  print '(8es12.4,/)',ave(1:ndim)
   do j=1,ndim
      print '(8es12.4)',(Sigma(j,i),i=1,ndim)
   end do

   temp = Sigma

   p=0

   call choldc(temp,p) !cholesky_decomposition

   L=0

   do i=1,ndim
      do j=i,ndim
         L(j,i) = temp(j,i)
      end do
      L(i,i) = p(i)
   end do

   print '(/,a)',' Choleksy L'
   do j=1,ndim
      print '(8es12.4)',(L(j,i),i=1,ndim)
   end do

   LT=transpose(L)

   print '(/,a)',' Choleksy LT'
   do j=1,ndim
      print '(8es12.4)',(LT(j,i),i=1,ndim)
   end do


   M = matmul(L,LT)

   print '(/,a)',' LLT'
   do j=1,ndim
      print '(8es12.4)',(M(i,j),i=1,ndim)
   end do

   do i=1,size(L,1)
      call ran_gauss(gauss(:)%vec(i))
   enddo

endif

call create_dist_with_correlation(L, muons, ave)

call compute_correlation(nmuons, muons, Sigma, ave)
print '(/,a)',' average and covariance matrix of new distribution'

print '(8es12.4)',ave(1:ndim)
print *
do j=1,ndim
   print '(8es12.4)',(Sigma(j,i),i=1,ndim)
end do

deallocate(p)
deallocate(Sigma)
deallocate(L,LT)
deallocate(M)
deallocate(temp)
deallocate(gauss)
deallocate(muons_ref)
deallocate(ave)

return
end subroutine match_covariance_distribution