	program sign_sim_mu
* 	Samvel Ter-Antonyan (Southern University, 10/03/2011)
C_____  Any change of simulation parameters: avL, sigmaL, stat0, sigma0, requres 
C_____  preliminary computation of "stot1" and "sigma1" parameters (see below) 

	PARAMETER (N=1000000, Mstr=85, Mdom=60)
	DOUBLE PRECISION v,u,W(99),s,sg,d,dg,stot0,x(99),y(99),xl(99)
	DOUBLE PRECISION AN,DOM,q,P,step,zmax,zmin,stot1,Ntot
	DOUBLE PRECISION sq,dq,epsilon,h0,h1,h2,dh0,dh1,dh2,r0,r1
	DOUBLE PRECISION sigma0,sigma1,sq0,dq0,sq1,dq1,x0,x1,x2
	DOUBLE PRECISION Mbgr0(200), Mbgr1(200), Mbgr2(200),vmin
	DOUBLE PRECISION Bank0(100,100),Bank1(100,100),stepv
	DOUBLE PRECISION StrA1(6,3),StrB1(6,3),StrA2(6,3),StrB2(6,3)
	DOUBLE PRECISION StrAB(6,3)
	INTEGER M(8),K0(320),K1(320),KG(320)
	REAL VECT(10000),kmut,sx(6),sy(6),dx(6),dy(6),rho(6)
	ND=Mdom*Mstr
	Ntot=real(N)*real(ND)
	zmin=-7.D0 ; zmax=7.D0 ; nbin=140 ; step=(zmax-zmin)/nbin
	vmin=100.D0 ; stepv=10.D0 ; kmut=0.00000000
	do i=1,6 ; do j=1,3                           !    Initializations of correlation arrays
 	  StrA1(i,j)=0.D0 ; StrB1(i,j)=0.D0
	  StrA2(i,j)=0.D0 ; StrB2(i,j)=0.D0 ; StrAB(i,j)=0.D0
	enddo; enddo
	do i=1,8 ; M(i)=0 ; enddo
	do i=1,nbin ; K0(i)=0 ; K1(i)=0 ; KG(i)=0 ; enddo
	do i=1,200;Mbgr0(i)=0.D0;Mbgr1(i)=0.D0;Mbgr2(i)=0.D0;enddo
C______ Simulation parameters
        Amu=4700.0 ; Dmu=sqrt(Amu)                    !  CORSIKA+GEANT for ICECUBE mu-hits
	avL=2.5 ; sigmL=0.7
C	avL=2.5 ; sigmL=0.87                          !  for UPPSALA presentation data

	stot0= 270.D0 !  input value
	sigma0=25.D0  !  input value
	stot1= 284.34381D0 !  is taken from previously simulated file "1c" for N=10^6
	sigma1=30.60764D0  !  is taken from previously simulated file "1c" for N=10^6

C        stot1= 286.39169D0 !  computed for UPPSALA fit and N=10^6
C        sigma1=30.07573D0  !  computed for UPPSALA fit and N=10^6

	AN=real(N) ; DOM=real(ND) ; sqd=dsqrt(DOM)

	sq0=0.D0 ; dq0=0.D0 ; sq1=0.D0 ; dq1=0.D0
	h0 =0.D0 ; h1 =0.D0 ; h2 =0.D0
	dh0=0.D0 ; dh1=0.D0 ; dh2=0.D0
	DO j=1,N                ! --------> Main event (timebase is 1 second) simulation loop
	  r0=0.D0 ; r1=0.D0
	  do ii=1,Mstr
	    do jj=1,Mdom
	      Bank0(ii,jj)=0.D0 ; Bank1(ii,jj)=0.D0
	    enddo
	  enddo
C___________________________________
	  do kstr=1,Mstr             ! ---- Filling all DOMs with Gausian noise
	    CALL RNORML(VECT,Mdom)   ! CERNLIB GAUSSIAN (avg=O,sigma=1) RANDOM NUMBER GENERATOR 
	    do kdom=1,Mdom
	      u=VECT(kdom)*sigma0+stot0
	      sg=sg+u ; dg=dg+u**2
	      Bank0(kstr,kdom)=u
	    enddo
	  enddo                      ! ----------------------------------------
 
	  CALL RNORML(vv,1) ; vv=vv*Dmu+Amu  !   random muon number
	  NPU=nint(vv)                       !   muons for 1 sec
	  
	  CALL RNORML(VECT,NPU)           
          do ii=1,NPU                !---> simulation of muon hits for DOM with index: jstr,id
	    v1=VECT(ii)*sigmL+avL    ! random Gaussian log(number) of mu-hits
            v1=exp(v1)
	    nv=nint(v1)
	    if(nv.eq.0) go to 8
	    if(nv.gt.ND) nv=ND
	    kmut=kmut+1.0
	    if(kmut.le.1000000.) write(70,*)nv
	    id=0                       
	    jstr=Mstr*rand()+1       ! random string
	    do kd=1,nv
	      id=id+1
	      if(id.gt.Mdom) then    ! immitation of depth dependence
	        id=1
	        jstr=Mstr*rand()+1
	      end if
	      Bank1(jstr,id)=Bank1(jstr,id)+1.D0
	    end do 
8 	  end do                     !-----------------------------------------------------

	  do kstr=1,Mstr             ! -- Data analysis --
	   do kdom=1,Mdom
	     x0=Bank0(kstr,kdom); h0=h0+x0 ; dh0=dh0+x0**2 !  No muons
	     x1=Bank1(kstr,kdom); h1=h1+x1 ; dh1=dh1+x1**2 !  Only Muon hits
	     x2=x0+x1 ; h2=h2+x2 ; dh2=dh2+x2**2           !  All hits
C___________ Correlation arrays______________________________
	     if(kstr.eq.10.or.kstr.eq.20.or.kstr.eq.30) then    ! (10,20,30 random selected 3 strings)
	       ks=kstr/10
	       if(kdom.eq.1. or.kdom.eq.11.or.kdom.eq.21.or.
     &            kdom.eq.31.or.kdom.eq.41.or.kdom.eq.51) then  ! these 6 numbers randomly selected
	          kd=(kdom-1)/10+1
	          y0=Bank0(kstr,kdom+1)
	          y1=Bank1(kstr,kdom+1)
	          y2=y0+y1
                  StrA1(kd,ks)=StrA1(kd,ks)+x2
                  StrB1(kd,ks)=StrB1(kd,ks)+y2
                  StrA2(kd,ks)=StrA2(kd,ks)+x2**2
                  StrB2(kd,ks)=StrB2(kd,ks)+y2**2
	          StrAB(kd,ks)=StrAB(kd,ks)+x2*y2
	       end if
	     end if
C_____________________________________________________________
             nb0=int((x0-vmin)/stepv+1.D0)
             nb1=int(x1+1.D0)
             nb2=int((x2-vmin)/stepv+1.D0)
             if(nb0.ge.1.and.nb0.le.nbin) Mbgr0(nb0)=Mbgr0(nb0)+1.D0
             if(nb1.ge.1.and.nb1.le.nbin) Mbgr1(nb1)=Mbgr1(nb1)+1.D0
             if(nb2.ge.1.and.nb2.le.nbin) Mbgr2(nb2)=Mbgr2(nb2)+1.D0
	     r0=r0+(x0-stot0) ; r1=r1+(x2-stot1)
	   enddo
	  enddo
	  r0=r0/DOM ; r1=r1/DOM
	  s0=r0/(sigma0/sqd) ; s1=r1/(sigma1/sqd)
	  nb0=int((s0-zmin)/step+1.0) ; nb1=int((s1-zmin)/step+1.0)
	  if(nb0.gt.1.and.nb0.le.nbin) K0(nb0)=K0(nb0)+1
	  if(nb1.gt.1.and.nb1.le.nbin) K1(nb1)=K1(nb1)+1
	  sq0=sq0+s0 ; dq0=dq0+s0**2 ; sq1=sq1+s1 ; dq1=dq1+s1**2
	  if(j.le.500000) write(71,88)s0,s1 
88 	  format(2f8.4)
	END DO !  -----------------------------------------------------


C-------OUTPUT---------------------------------------------------------
	OPEN(90,file='0c')
	OPEN(91,file='1c')
	OPEN(92,file='2c')
	open(93,file='3mc')
	v0=vmin+stepv/2.0
	v2=vmin+stepv/2.0
	do i=1,nbin
	  write(90,91)v0,Mbgr0(i)/real(ND)
	  write(92,91)v2,Mbgr2(i)/real(ND)
	  v0=v0+stepv
	  v2=v2+stepv
	end do
	
	do i=1,nbin
	  write(91,91)real(i),Mbgr1(i)/real(ND)
	end do

	h0=h0/Ntot ; dh0=dsqrt(dh0/Ntot-h0**2)
	h1=h1/Ntot ; dh1=dsqrt(dh1/Ntot-h1**2)
	h2=h2/Ntot ; dh2=dsqrt(dh2/Ntot-h2**2)

	write(91,*)' '
	write(91,*)N
	write(91,*)h0,dh0
	write(91,*)h1,dh1
	write(91,*)h2,dh2

	write(*,*)h0,dh0
	write(*,*)h1,dh1
	write(*,*)h2,dh2

	z=zmin+step/2.0
	do i=1,nbin
	  write(93,20)z,K0(i),K1(i)    !    Significance distributions 
	  z=z+step
	end do
	write(93,*)

	sq0=sq0/AN ; dq0=dsqrt(dq0/AN-sq0**2)
	sq1=sq1/AN ; dq1=dsqrt(dq1/AN-sq1**2)
	write(93,*)'NTot',N
	write(93,*)'0mu:',sq0,dq0
	write(93,*)'All:',sq1,dq1
	write(*,*) 'Significance'
	write(*,*) '0mu:',sq0,dq0
	write(*,*) 'All:',sq1,dq1
	
	do is=1,3
	  do id=1,6
	    sx(id)=StrA1(id,is)/AN
	    sy(id)=StrB1(id,is)/AN
	    dx(id)=sqrt(StrA2(id,is)/AN-sx(id)**2)
	    dy(id)=sqrt(StrB2(id,is)/AN-sy(id)**2)
	    sxy=StrAB(id,is)/AN
	    rho(id)=(sxy-sx(id)*sy(id))/(dx(id)*dy(id)) !  correlations for 6 depths
	  end do
	  write(50,*)is
	  write(50,50)rho                   
	  write(51,*)is
	  write(51,50)sx
	  write(51,*)''
	  write(51,50)dx
	  write(52,*)is
	  write(52,50)sy
	  write(52,*)''
	  write(52,50)dy
	end do

	stop 
20	format(f10.3,3i12)
50	format(6f12.4) 
91	format(f12.1,3e12.5) 
93	format(3i12)
	end
