!plot using g-2/documentation/EfieldSystematic/start.gnu
program generate_electron_dist
  use precision_def
  implicit none
  integer i,j
  integer number_per_ybin(0:101), iy_bin, Nbins/100/
  integer number_per_tbin(0:401,0:401)  
  integer tbin

  real(rp) angle(0:401)
  real(rp) theta, energy,phi
  real(rp) y
  real(rp) pi/3.141592/
  real(rp) ef(0:101)

number_per_ybin(:)=0
number_per_tbin(:,:)=0

do i=1,1000000
 call electronDist(energy,theta,phi)
 write(111,'(3es12.4)')energy,theta,phi

  y= energy/52.83e6

  iy_bin = y * Nbins +.5
  ef(iy_bin) = y
  number_per_ybin(iy_bin)= number_per_ybin(iy_bin)+1  !histogram energy
!histogram theta for the energy bin
  tbin = theta * Nbins/pi + .5
  number_per_tbin(iy_bin, tbin) = number_per_tbin(iy_bin, tbin)+ 1
  angle(tbin) = theta
end do
do i=1,nbins
write(112,'(2i10)')i,number_per_ybin(i)
end do
do i=1,nbins
  do j=1, nbins
    write(113,'(2i10,1x,2es12.4,i10)')i,j,ef(i),angle(j),number_per_tbin(i,j)
  end do
end do


end