program Pair_DOMs
*       Samvel Ter-Antonyan (Southern University, 10/03/2011)
	DOUBLE PRECISION s1,s2,s12,d1,d2,r1,r2,w
        CALL POISET(20.)                   ! Initialization of CERNLIB Poisson random number generator
	N=1000000 ; w=real(N)              ! Total number of events (seconds) 
	av=270.D0 ;                        ! Average of Poisson noise
	mu=0			           ! Current average Poisson number of muons (muon hits) 
	DO k=1,150                         ! Loop for current average number of muons (muon hits) per DOM	             
	  write(*,*)k
	  s1=0.D0 ; d1=0.D0 ; s12=0.D0     ! Initializations of summators 
	  s2=0.D0 ; d2=0.D0                ! Initializations of summators
	  do i=1,N                         ! Loop of events (seconds)   
	    
	    CALL POISSN(av,Noise1,ierr)    ! Poisson noise of the fisrt  DOM   
	    CALL POISSN(av,Noise2,ierr)    ! Poisson noise of the second DOM
	    r1=Noise1 ; r2=Noise2

	    CALL POISSN(real(mu),Nmu,ierr) ! Simulated Poisson number of muons, Nmu (or muon hits)
	    r1=r1+Nmu                      ! Here is a correlation
	    r2=r2+Nmu                      ! Here is a correlation
	    s1=s1+r1
	    d1=d1+r1**2
	    s2=s2+r2
	    d2=d2+r2**2
	    s12=s12+r1*r2
	  end do
C_________Statistical analysis______________
	  avr1=s1/w                         ! average noise
	  rms1=dsqrt(d1/w-avr1**2)          ! RMS
	  avr2=s2/w
	  rms2=dsqrt(d2/w-avr2**2)
	  rho=(s12/w-avr1*avr2)/(rms1*rms2) ! Correlation coefficient
	  drho=sqrt((1.0-rho**2)/(w-2.D0))  ! Error of correlation coefficient
	  write(45,10)mu,rho,drho
	  mu=mu+2                           ! Increment for current average muon (hit) number  
	END DO
	stop
10	format(i7,5e12.4)
	end