subroutine compute_correlation(nmuons, muons, Sigma, ave )

use muon_mod
implicit none

type (params), allocatable ::  param(:)
type (muon_struct), allocatable :: muons(:)

real(rp), allocatable :: Sigma(:,:) 
real(rp), allocatable :: ave(:)
integer nmuons, i,j,k
integer ndim

ndim = size(Sigma,1)

if(.not. allocated(param))then
   allocate(param(1:nmuons))
endif
do i=1,nmuons
   param(i)%vec(1:6)=muons(i)%coord%vec(1:6)
   param(i)%vec(7)=muons(i)%coord%spin(1)
   param(i)%vec(8)=muons(i)%coord%spin(3)
end do

do j=1,ndim
   ave(j) = sum(param(:)%vec(j))/nmuons
   do i=1,ndim 
      Sigma(i,j) = sum( param(:)%vec(j) * param(:)%vec(i))/nmuons
   end do
end do

do i=1,ndim
   do j=1,ndim
      Sigma(i,j) = Sigma(i,j) - ave(j)*ave(i)
!      Sigma(i,j) = Sigma(i,j)/sqrt(Sigma(i,i)*Sigma(j,j))
   end do
end do


return
end subroutine compute_correlation



subroutine create_dist_with_correlation(L, muons, ave)
use muon_mod
implicit none
type (muon_struct), allocatable :: muons(:)
type (params), allocatable :: param(:), gauss(:)
real(rp), allocatable :: L(:,:), ave(:)
real(rp) random/6/
integer i,ndim,nmuon
integer k

nmuon = size(muons)
ndim=size(L,1)

allocate(gauss(size(muons)))
allocate(param(size(muons)))

do i=1,ndim
   call ran_gauss(gauss(:)%vec(i))
enddo
do i=1,nmuon
   param(i)%vec(:) = matmul(L,gauss(i)%vec(:) )
end do
do i=1,ndim
   if(i<=6)muons(:)%coord%vec(i)=param(:)%vec(i) + ave(i)
   if(i == 7)muons(:)%coord%spin(1)=param(:)%vec(i) + ave(i)
   if(i == 8)muons(:)%coord%spin(3)=param(:)%vec(i) + ave(i)
end do
muons(:)%coord%spin(2) = sqrt(abs(1.- muons(:)%coord%spin(1)**2- muons(:)%coord%spin(3)**2))

deallocate(gauss, param)
return
end subroutine create_dist_with_correlation



subroutine get_number_events(unit,  tot)

use muon_mod
implicit none
integer unit
integer i, tot,j, ios
integer nread
integer initx, inity, initz
character*290 line
character(72) plotting_script(4)/'hist_initial_dist.gnu','hist_after_inflector.gnu','hist_start_tracking.gnu','hist_out.gnu'/
logical include_spin/.false./
logical include_spin_valetov/.false./

type (muon_struct), allocatable :: muons(:)
allocate(muons(1:1))



tot=1
nread=0
do while(.true.)
  read(unit,'(a)', end=99) line
  if((index(line(1:3),'!')/= 0 .or. index(line(1:3),'#')/= 0) .and. (index(line,'s_x')/=0))include_spin=.true.
  if((index(line(1:3),'!')/= 0 .or. index(line(1:3),'#')/= 0) .and. (index(line,'sx')/=0))include_spin_valetov=.true.  
  if(index(line(1:3),'!')/= 0 .or. index(line(1:3),'#')/= 0)cycle
  if(line(1:10)== '          ')cycle

  i=tot

  i=1
  if(include_spin_valetov)then
          read(line,*, end=99) j, muons(i)%Jx, muons(i)%Jy, muons(i)%coord%vec(1), &
           muons(i)%coord%vec(2),  muons(i)%coord%vec(3), muons(i)%coord%vec(4), muons(i)%coord%vec(5), muons(i)%coord%vec(6),&
           muons(i)%coord%t, muons(i)%coord%s, muons(i)%coord%spin, initx, inity, initz
    elseif(include_spin)then
          read(line,'(i10,1x,13es16.8)') j, muons(i)%Jx, muons(i)%Jy, muons(i)%coord%vec(1), &
           muons(i)%coord%vec(2),  muons(i)%coord%vec(3), muons(i)%coord%vec(4), muons(i)%coord%vec(5), muons(i)%coord%vec(6),&
           muons(i)%coord%t, muons(i)%coord%s, muons(i)%coord%spin
     else 
           read(line,'(i10,1x,9es16.8)') i, muons(i)%Jx, muons(i)%Jy, muons(i)%coord%vec(1), &
           muons(i)%coord%vec(2),  muons(i)%coord%vec(3), muons(i)%coord%vec(4), muons(i)%coord%vec(5), muons(i)%coord%vec(6),&
           muons(i)%coord%t
  endif

  tot = tot+1
  
end do
99 print '(a,i10,a)',' Skip ', nread,' muons'
     print '(/,a,1x,i10,1x,a,1x,i10)','Phase space vector for ',tot,' particles read from unit ',unit
return
end subroutine get_number_events