fort/file 20 kdev51.ps
meta 20 -111
ve/cre ibur(12) 
ve/cre irun(12)
ve/cre ispi(12) 
ve/cre ispf(12) 
ve/cre inoi(12)
ve/cre n100(12)
ve/cre irms(12)
ve/cre irms1(12)
ve/cre icc(1)
ve/cre cont2(231,500)
ve/cre ccont2(231,500)
ve/cre cont3(231,50)
ve/cre ccont3(231,50)
ve/cre kstore(500000) i
ve/cre hstore(500000) r
ve/cre scal(1)
ve/cre imx5(1)
ve/cre imn5(1)
ve/cre vhskew(1)
ve/cre vhkur(1)
ve/cre ifreq(1)
ve/cre vhskew1(1)
ve/cre vhkur1(1)
ve/cre vfrms(1)
ve/cre vcrem(1)
do k =1,8
vec/inp icc [k]

application comis quit
c234567
      subroutine spillstat4
      common/pawc/hmemor(1000000)
      dimension cont(8000),aamp(2000),ccont(1000)
c      dimension cont2(231,50),ccont2(231,50)
      dimension rs(2000),nhsta(1)
      character*6 aa
      character*60 sfile
c     double precision smean5,rms5,rms52
      vector ibur,irun,ispi,inoi,ispf,n100,irms,irms1
      vector icc,scal,imx5,imn5,vhskew,vhkur,ifreq
      vector vhskew1,vhkur1,vfrms,vcrem
      vector cont2
      vector ccont2    
      vector cont3
      vector ccont3    
      vector kstore,hstore

c      call hlimit(2000000) 
      call hbook1(30110,' spill 5ms ',1000,1000.,6000.,0.) ! full range
      call hbook1(31110,' hits for spill 5ms ',1000,1000.,6000.,0.) ! full range
      call hminim(30110,1.0)
      call hminim(31110,1.0)
      call hbook1(32110,' hits / spill 5ms ',1000,1000.,6000.,0.) ! full range
      call hbook1(33110,' spill / hits 5ms ',1000,1000.,6000.,0.) ! full range
      call hbook1(110,' spill 5ms ',1000,1000.,6000.,0.) ! full range
      call hbook1(1110,' hits forspill 5ms ',1000,1000.,6000.,0.) ! full range
      call hmaxim(110,1000.)
      call hmaxim(33110,2.5)

      call hbook1(101,' spill 1ms ',4000,1500.,5500.,0.)
      call hbook1(102,' spill 5ms ',800,1500.,5500.,0.)
      call hbook1(103,' spill 25ms ',160,1500.,5500.,0.)
      call hbook1(104,' spill .5 ms ',1000,1500.,2000.,0.)

      call hbook1(200,' mean 1 ms ',50,50.,150.,0.)
      call hbook1(201,' sd/sqrt(mean) 1 ms ',60,1.,4.,0.)
      call hbook1(202,' skewness 1 ms ',50,-1.,1.,0.)
      call hbook1(203,' kurtosis  1 ms ',40,1.,5.,0.)
    
      call hbook1(204,' lf noise ',25,0.,0.5,0.)
      call hbook1(205,' 100 Hz noise ',25,0.,0.25,0.)
      call hbook1(206,' total  noise ',30,0.,0.6,0.)
      call hbook1(207,' high spill % ',25,0.,25.,0.)
      call hbook1(307,' high folded spill % ',50,0.,50.,0.)
      

 
      call hbook2(301,' sd/sqrt(mean) vs skewness ',
     1 60,1.,4.,50,-1.,1.,0.)
      call hbook2(302,' sd/sqrt(mean) vs kurtosis ',
     1 60,1.,4.,40,1.,5.,0.)
      call hbook2(303,' skewness  vs kurtosis ',
     1 50,-1.,1.,40,1.,5.,0.)

      call hbook2(304,' sd/sqrt(mean) vs lf noise ',
     1 60,1.,4.,25,0.,0.5,0.)
      call hbook2(305,' sd/sqrt(mean) vs 100 Hz',
     1 60,1.,4.,25,0.,0.25,0.)
      call hbook2(306,' sd/sqrt(mean) vs total noise',
     1 60,1.,4.,30,0.,0.6,0.)

      call hbook1(800,' spill - folded ',231,0.,924.,0.) ! for sps time det

      call hbook1(8000,' freq ',250,0.,250.,0.)
      call hbook1(9000,' RMS ',50,0.,1000.,0.)
      call hbook1(9100,' First Difference 1 ms ',50,-100.,150.,0.)
      call hbook1(9200,' correlelogram',2000,0.,1000.,0.)
      call hbook2(10000,' corrected folded spill, spill ',
     1 231,0.,924.,500,1100.,6100.,0.)
      call hbook1(11000,' Corrected folded time',
     1 231,0.,924.,0.)
      call hbook1(12000,' Uncorrected folded time',
     1 231,0.,924.,0.)
      call hbook1(12100,' first diff.Uncorrected folded time',
     1 231,0.,924.,0.)
      call hbook1(12300,' 2ndt diff.Uncorrected folded time',
     1 231,0.,924.,0.)
      call hbook1(12200,' fisrt diff.Uncorrected folded time distr.',
     1 150,-500.,1000.,0.)

      call hbook1(7000,' hits ',125,0.,250.,0.)
      call hbook1(7100,' hits 2.0 - 5.0 ',125,0.,250.,0.)
      call hbook2(20000,'  spill, hits ',
     1 500,1100.,6100.,50,0.,250.,0.)
      call hbook2(21000,' folded spill, hits ',
     1 231,0.,924.,50,0.,250.,0.)

       if(icc(1).eq.1)sfile = 'elis_2022-2567-0336-ian02.txt'
       if(icc(1).eq.2)sfile = 'elis_2022-2066-1129-ian02.txt'
       if(icc(1).eq.3)sfile = 'elis_2022-2288-1360-ian02.txt'
       if(icc(1).eq.4)sfile = 'elis_2022-2567-1146-ian02.txt'
       if(icc(1).eq.5)sfile = 'elis_2023-3165-0035-ian02.txt'
       if(icc(1).eq.6)sfile = 'elis_2023-3165-0049-ian02.txt'
       if(icc(1).eq.7)sfile = 'elis_2023-3165-0056-ian02.txt'
       if(icc(1).eq.8)sfile = 'elis_2023-3229-0170-ian02.txt'


    
       open(unit=19,file=sfile,status='unknown') 


*      if(icc(1).eq.2)
*    1 open(unit=19,file='elis-12567-0337.txt',status='unknown') ! golden file 
*      if(icc(1).eq.3)
*    1 open(unit=19,file='elis-12567-0338.txt',status='unknown') ! golden file 
*      if(icc(1).eq.4)
*    1 open(unit=19,file='elis-12567-0339.txt',status='unknown') ! golden file 
*      if(icc(1).eq.5)
*    1 open(unit=19,file='elis-12567-0340.txt',status='unknown') ! golden file 
*      if(icc(1).eq.6)
*    1 open(unit=19,file='elis-12567-1146.txt',status='unknown') ! golden file 
*      if(icc(1).eq.7)
*    1 open(unit=19,file='elis-2465-0170.txt',status='unknown') ! dataplot31.kumac
*      if(icc(1).eq.8)
*    1 open(unit=19,file='elis-2566-1230.txt',status='unknown') ! dataplot32.kumac
*      if(icc(1).eq.9)
*    1 open(unit=19,file='elis-2066-1129.txt',status='unknown') ! dataplot33.kumac
*     if(icc(1).eq.10)
*    1 open(unit=19,file='elis-12288-1360.txt',status='unknown') ! gold 22
*     if(icc(1).eq.11)
*    1 open(unit=19,file='elis-13019-0133.txt',status='unknown') !  new 23
*     if(icc(1).eq.12)
*    1 open(unit=19,file='elis-13046-0011.txt',status='unknown') ! new 23

      read(19,5510,END=999)aa,nrun,nbur
*5510 format(a6,2i8)
 5510 format(a6,2i9)
      write(*,*)' '
      write(*,*)' run burst ',nrun,nbur
      ii = icc(1)
      ibur(1) = nbur
      irun(1) = nrun

      nnn =0
      nn = 500000 !  max input - see arrays
      do k =1,nn
      read(19,5513,END=999)aa,kk,kid,kcode,kstamp,nhsta(1)
*     if(nhsta(1).gt.100.)go to 444 
      nnn = nnn+1
 5513 format(a6,2i9,2z10,i8)
      kstore(nnn) = kstamp !  for ffb determination
      hstore(nnn) = nhsta(1) !  store hits for folded plots
      fm = float(kstamp)/40000d0 ! ms
      call hfill(110,fm,0.,1.) ! full spill 5 ms
      call hfill(30110,fm,0.,1.) ! full spill 5 ms
      temp = nhsta(1)
      call hfill(31110,fm,0.,temp) ! full spill 5 ms hits sum
      call hfill(101,fm,0.,1.)
      call hfill(102,fm,0.,1.)
      call hfill(103,fm,0.,1.)
      call hfill(104,fm,0.,1.)

      temp = nhsta(1)
      call hfill(7000,temp,0.,1.) ! hits
      if(fm.gt.2000..and.fm.lt.5000.)then
      call hfill(7100,temp,0.,1.) ! hits
      endif
c     call hfill(20000,temp,fm,1.) ! hits vs spill
      call hfill(20000,fm,temp,1.) !  spill vs hits
c--------------------------------
       go to 444
c     folded data

      ffb = 923.9926  !  fold time approx
      tempx = mod(float(kstamp),ffb)  !  folded spill
      tempy = float(kstamp)/40000d0  ! ms spill
      xxx = (tempy -1500.)/1000.
      corr = 2.4*xxx**3
      ctempx = tempx - corr ! corrected folded spill
      if(ctempx.lt.0.)ctempx = ctempx +924.0
*      call hfill(21000,ctempx,temp,1.) ! corrected folded spill vs hits
c     note that above plot is reordered ( see below) 
c     but an approx sps time is used here. <<<<<<<<<<<<<<NOTE<<<<<<<<<
c      This 
c     avoids  storing the GTK hits as was done for kstamp.
c     The consequence is that the reordering is imperfect.<<<<<<<<<<<<

*     call hfill(10000,ctempx,tempy,1.)
*     call hfill(11000,ctempx,0.,1.)
*     if(tempy.gt.1200..and.tempy.lt.2200.)then
c     uncorrected folded
*     call hfill(12000,tempx,0.,1.)
*     endif

 444  continue
c--------------------------------
      enddo
 999  continue

      write(*,*)' count  ', nnn
      write(*,*)' '
      close (19)
c------------------------------------------------------------------

c     ratio hits/triggers
      c1 =hsum(31110)
      c1 =1./c1
      c2 = hsum(30110)
      c2 = 1./c2
      call hopera(31110,'/',30110,32110,c1,c2)
      call hopera(30110,'/',31110,33110,c2,c1)


c------------------------------------------------------------------
c      determine sps time  *****************************************************
      
c     goto 555
 800  continue 
      call hreset(800,' ')    
      ffs = 923.9910 ! mod 24/2/23
      ffb=  0
      cmin2 = 5000.
      do i=1,50 ! mod 24/2/23 
      ff = i
      fff = ffs+(ff-1.)*0.00005
      do k=1,nnn
      dtemp = mod(float(kstore(k)),fff) ! change  26/2/23
c     correct the time here -----------------------------
c     4th order correction - apply over full range
      tempy = kstore(k)/40000.d0 
      xxx = tempy - 1500.
      xxx = xxx/1000.
      corr = 2.4*xxx**3    ! xxx in sec 
      temp = dtemp
      temp=temp-corr
      if(temo.lt.0.)temp = temp +924. ! added here 3/6/23
C     if(k.lt.20)write(*,*)' ',k,i, kstore(k),temp 

      call hfill(800,temp,0.,1.)
      enddo  !k
      call hunpak(800,cont,' ',0)
c     write(*,*)' fff ',fff
c     write(*,*)' cont ',(cont(kk),kk= 1,231)
      cmin = 5000.
      do kk = 2,230
      if(cont(kk).lt.cmin)then
      cmin = cont(kk)
      kkmin = kk
      endif
      enddo
      call hix(800,kkmin,var) 
      write(*,*)' fff , cmin, kkmin ', fff,cmin,kkmin 
      if(cmin.lt.cmin2)then
      cmin2 = cmin
      ffb = fff
      endif
      call hreset(800,' ')
      enddo  !
      scal(1) = ffb - 923.
      write(*,*)' --------------------------'
      write(*,*)' '
      write(*,*)' optimum value ',ffb 
      write(*,*)' '
      write(*,*)' --------------------------'
********************** *******************************

 555  continue 




c     get folded plots using optimised sps time - ffb
      do k = 1,nnn 
      kstamp = kstore(k)
      tempx = mod(float(kstamp),ffb)  !  folded spill
      tempy = float(kstamp)/40000d0  ! ms spill
      xxx = (tempy -1500.)/1000.
      corr = 2.4*xxx**3
      ctempx = tempx - corr ! corrected folded spill
      if(ctempx.lt.0.)ctempx = ctempx +924.0
      call hfill(10000,ctempx,tempy,1.) ! fig 1
      call hfill(11000,ctempx,0.,1.) !  not used
      temp = hstore(k) !   optimised sps ffb used here for folded spill
      call hfill(21000,ctempx,temp,1.) ! corrected folded spill vs hits Fig 7
      if(tempy.gt.1200..and.tempy.lt.2200.)then ! standard
*     if(tempy.gt.2000..and.tempy.lt.5000.)then ! temp mod
c     uncorrected folded
      call hfill(12000,tempx,0.,1.) ! used - Fig 3
      endif
      enddo ! k loop







c------------------------------------------------------------------


c     go to 333  !  do not reorder <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

c     REORDER time plots
c     comis will not work with the following code
c     see spillstat5.f  for working code
c     cont2, ccont2 require too much memory.
c     but OK with a vector to compile.
c     conclude vectors OK, dimension in comis fortran fails.
c     and note that these differ from kstat44.kumac
c     due to different sps fold time.  Also
c     ordering may change since different minimum is selected.
c     
c     reorder folded spill plot..................................
      call hunpak(12000,cont,' ',0) ! projection
      call hunpak(10000,cont2,' ',0) ! 2D spill
      call hunpak(21000,cont3,' ',0) ! 2D hits
      cmin = 1000000.
      jmin = 1
      do j =1,231
      if(cont(j).lt.cmin)then ! temp remove
      cmin = cont(j)
      jmin = j
      endif
      enddo 
      do j = 1,231
      k =  j - jmin + 1
      if(k.lt.1)k = k +231
      ccont(k) = cont(j)
      do jj = 1,500
      ccont2(k,jj) =cont2(j,jj)
      enddo
      do jj =1,50
      ccont3(k,jj) = cont3(j,jj) 
      enddo
      enddo
      call hpak(12000,ccont)
      call hpak(10000,ccont2)
      call hpak(21000,ccont3)
c     end reorder
 333  continue


c     examine 2D spill / folded spill rordered
      fmaxcont = 0
      fmean = 0.
      frms = 0.
      crem = 0.
      do j = 1,231
      do jj = 1,500
      if(ccont2(j,jj).gt.fmaxcont)fmaxcont = ccont2(j,jj)
      if(ccont2(j,jj).gt.0.)then  !4
      fmean = fmean +ccont2(j,jj)
      frms = frms + ccont2(j,jj)**2 
      fnum = fnum + 1.
c     remove bad regions <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      jk = 0
      jjk =0
      csave = 0
*     if(ccont2(j,jj).gt.20.)then !3
      if(ccont2(j,jj).gt.10.)then !3  240623
      csave = ccont2(j,jj) ! triggers removed
      crem = crem +csave ! sum removed triggers
      ccont2(j,jj) = 0. ! <<<<<<<<<<<<<<<< remove 2d
      jk = j
      jjk = jj
      endif !3
      if(jk.ne.0.and.jjk.ne.0)then !2
c     bad area
      call hijxy(10000,jk,jjk,x,y) ! edge of bad region
c     correct 12000 (ccont) if in spill range
      if(y.gt.1200.0.and.y.lt.2200.)then ! 1
      ccont(jk) = ccont(jk) - csave    ! remove 1d in spill range
      endif ! 1
      endif ! 2
      endif ! 4
      enddo 
      enddo
c**********************************REPLACE******************
      call hpak(10000,ccont2)
      call hpak(12000,ccont)
c     bad regions removed 
c**********************************REPLACE******************

      frms = sqrt(frms/fnum)
      fmean = fmean/fnum
      crem = crem/float(nnn)*100.
      vcrem(1) = crem
      write(*,*)' ------------------------'
      write(*,*)' ------------------------'
      write(*,*)' ccont2 mean rms max ',fmean,frms, fmaxcont
      write(*,*)' % spill removed ',crem
      write(*,*)' ------------------------'
      write(*,*)' ------------------------'
     
      
c     ---------------------------------------
c     folded uncorrected spill   statistics
      call hunpak(12000,cont,' ',0)
      tot = 0.
      f500 = 0.
      smean = 0
      fj = 0.
      do j = 1,231
      if(cont(j).gt.200.)then
      smean = smean + cont(j)
      fj = fj  +1.
      endif
      enddo
      smean = smean/fj
      rms = 0
      do j = 1,231
      if(cont(j).gt.200.) rms = rms +(cont(j) -smean)**2
      enddo
      rms = sqrt(rms/fj)
      itemp =rms
      irms(1) = itemp
      write(*,*)' ---------------------------'
      write(*,*)' RMS folded spill ',ii,rms
      write(*,*)' ---------------------------'
      do j =1,231
      tot = tot + cont(j)
      if(cont(j).gt.500.)f500 = f500 +cont(j)
      enddo
      f500 = f500/tot  !fraction gt 500
      f500 = f500*100.
      itemp = f500
      ispf(1) = itemp
      call hfill(307,f500,0.,1.)
c--------------------------------------------------------
c     SPILL plots
      call hunpak(102,cont,' ',0) ! 5 ms interval 1.5 - 5.5
      tot = 0.
      t600 =0.
      do jj =1,800
      call hfill(9000,cont(jj),0.,1.)
      tot = tot + cont(jj)
      if(cont(jj).gt.600.)t600 = t600 + cont(jj)
      enddo
      write(*,*)' '
      f600 = t600/tot ! fraction entries gt 600
      itemp = f600*100
      ispi(1) = itemp
c     1.5 - 5.5 sec
      call hunpak(101,cont,' ',0) ! 1 ms interval

c     first difference plots
      do i = 1,3999
      temp = cont(i+1) -cont(i) ! first diference of spill
      call hfill(9100,temp,0.,1.)
      enddo




c     spill statistics for 1 ms
      fnn = 4000
      smean = 0.
      do i = 1,fnn
      smean = smean +cont(i)
      enddo
      smean = smean/fnn
      sd = 0.
      skew = 0.
      skurt = 0.
      do i = 1,fnn
      sd = sd + (cont(i) - smean)**2
      skew = skew + (cont(i) - smean)**3
      skurt = skurt + (cont(i) - smean)**4
      enddo

      sd =sqrt(sd/fnn)
      isd = sd
      irms1(1) = isd  !  for display 1ms
      skew = skew/sd**3/fnn
      skurt = skurt/sd**4/fnn
c     try excess RMS <<<<<<<<<<<<<<<<<<
      sd = sqrt(sd**2-n)  !  mod
      sdn = sd/sqrt(smean)

      write(*,*)' mean,sd/sqrt(mean),skewness,kurtosis ',
     1 ii, smean,sdn,skew,skurt    
      write(*,*)' '
      write(*,*)' %fraction gt 600.',ii,f600*100.
      write(*,*)' %folded fraction gt 500.',ii,f500
      f600 = f600*100. !  to %



      call hfill(200,smean,0.,1.)
      call hfill(201,sdn,0.,1.)
      call hfill(202,skew,0.,1.)
      call hfill(203,skurt,0.,1.)
      call hfill(207,f600,0.,1.)     

      call hfill(301,sdn,skew,1.)
      call hfill(302,sdn,skurt,1.)
      call hfill(303,skew,skurt,1.)
c-------------------------------------------


*     periodogram    
 
      ssm = 0.
      do i =1,1000
      ssm = ssm + cont(i)
      enddo
      ssm = ssm/1000.
*     ssm = 0. !    <<<<<<<<< note <<<<<<<<<<<<<
      do i = 1,1000
      cont(i) = cont(i) - ssm
      enddo
      pi2 = 6.283185
      pi = 3.141592654
c     form periodogram
      do iw = 0,999 ! freq
      xs = 0.
      xc =0.
      do i = 1,1000
      t = i
      t = t/1000.  ! ms
      fq = iw
      xs = xs + cont(i)*sin(pi2*fq*t)
      xc = xc + cont(i)*cos(pi2*fq*t)
      enddo
      aamp(iw+1) = (xs**2 + xc**2)
      aamp(iw+1) = sqrt(aamp(iw+1))
      temp =iw
      call hfill(8000,temp,0.,aamp(iw+1))
      enddo
c-------------------------------------------------
c     correlogram

      call hunpak(104,cont,' ',0)
      yb = 0
      ny = 0. 
      il =1
      ih = 1000
      do i = il,ih
      yb = yb + cont(i)
      ny = ny +1
      enddo

      yb = yb/float(ny)
      y2 = 0.
      do i = il,ih
      y2 = y2 + (cont(i)-yb)**2
      enddo



      write(*,*)' '
      do is = 1,ny-1
      rs(is) = 0.
      do it = 1,ny -is
      rs(is) = rs(is) + (cont(it)-yb)*(cont(it+is)-yb)/y2
*     write(*,*)' is it it+is ',is, it , is +it
      enddo
*     write(*,*)' spill is rs ', is,'  ',rs(is)
      temp = is 
      temp = temp*0.5
      temp1= rs(is)*float(ny)/float(ny-is) 
      call hfill(9200,temp,0.,temp1)
      enddo 
*     write(*,*)' '

*     periodogram   end 
c------------------------------------------------
c     noise levels
c     background ~ 400.
      sn1 = 0.
      do j = 1,40
      sn1 = sn1 + aamp(j) - 400.
      enddo
      sn2 = 0.
      do j = 95,105
      sn2 = sn2 + aamp(j) - 400.
      enddo
c     sn2 = abs(sn2)  ! avoid small negative numbers
      sn3 = 0.
      do j = 1,155
      sn3 = sn3 + aamp(j) !  total signal
      enddo
      sn4 = 0.
      do j = 1,155
      sn4 =  sn4 + aamp(j) -400.  ! total noise 
      enddo

c     max amplitude
      amax = 0
      do j = 1,500
      if(aamp(j).gt.amax)then
      jmax = j
      amax = aamp(j)
      endif 
      enddo       
     
c     for consistency - less sensitive to 400 subtraction
c     but can omit 50 Hz signal

c     sn4 = sn1 + sn2        

      noise = sn4*100./sn3 ! total noise
      inoi(1) = noise
      noise = sn2/sn3*100. !  100 Hz
      n100(1)  = noise
c      replace by max amplitude
      itemp = amax
      n100(1) = itemp 
      ifreq(1) = jmax -1 
      write(*,*)' sn2-100,sn4-tot-noise,sn3-tot-sig',
     1  sn2,sn4,sn3
      write(*,*)' noise:  low ,100 Hz, total ',sn1,sn2,sn3        
      write(*,*)' noise fractions low, 100 Hz, total noise ',
     1  sn1/sn3, sn2/sn3, sn4/sn3

      call hfill(204,sn1/sn3,0.,1.)
      call hfill(205,sn2/sn3,0.,1.)
      call hfill(206,sn4/sn3,0.,1.)

      call hfill(304,sdn,sn1/sn3,1.)
      call hfill(305,sdn,sn2/sn3,1.)
      call hfill(306,sdn,sn4/sn3,1.)

c234567
c     maximum in 5ms spill plot
c     + mean, rms  range limited
      call hunpak(110,cont,' ',0)
      smax = 0.
      smin = 1000.
      smean5 = 0.
      rms5 = 0.
      fn5 =0.
      do j = 1,1000
      if(cont(j).gt.smax)smax = cont(j)
      if(j.gt.200.and.j.lt.800)then
      if(cont(j).lt.smin)smin = cont(j)
      if(abs(cont(j) - 500.).lt.250.)then
      smean5 = smean5 + cont(j)
      rms5 = rms5 + cont(j)**2
      fn5 = fn5 + 1. 
      endif
      endif
      enddo
      imx5(1) = smax
      imn5(1) = smin
      smean5 = smean5/fn5
      rms5 =rms5/fn5
      rms5  = sqrt(rms5 -smean5**2)

c     check
      rms52 = 0.
      fn5 =0.
      do j =1,1000
      if(j.gt.200.and.j.lt.800)then
      if(abs(cont(j) - 500.).lt.250.)then
      rms52 = rms52 + (cont(j)-smean5)**2
      fn5 = fn5 + 1. 
      endif
      endif
      enddo
      rms52 = sqrt(rms52/fn5)

      write(*,*)' '
      write(*,*)' MIN/MAX       ',ii,smin,smax 
      write(*,*)' 5 ms mean rms ',ii,smean5,rms5,rms52
      write(*,*)' '

c     stats for hits

      call hunpak(7000,cont,' ',0)
      hmean = 0.
      fn = 0.
      do j =1,125
      v = float(j)*2. -1. ! bin center
      hmean = hmean + cont(j)*v
      fn = fn +cont(j)
      enddo
      hmean = hmean/fn
      hrms = 0.
      hskew = 0.
      hkur = 0.
      do j =1,125
      v = float(j)*2. -1. ! bin center
      hrms = hrms + cont(j)*(v - hmean)**2
      hskew = hskew + cont(j)*(v - hmean)**3
      hkur = hkur + cont(j)*(v - hmean)**4
      enddo
      hrms = sqrt(hrms/fn)
      hskew = hskew/fn/hrms**3
      hkur = hkur/fn/hrms**4
      write(*,*)' '
      write(*,*)' HITS , mean rms  skewness kurtosis',
     1 ii,hmean,hrms,hskew,hkur
      write(*,*)' '
      vhskew(1) = hskew
      vhkur(1) = hkur


      call hunpak(7100,cont,' ',0)
      hmean = 0.
      fn = 0.
      do j =1,125
      v = float(j)*2. -1. ! bin center
      hmean = hmean + cont(j)*v
      fn = fn +cont(j)
      enddo
      hmean = hmean/fn
      hrms = 0.
      hskew = 0.
      hkur = 0.
      do j =1,125
      v = float(j)*2. -1. ! bin center
      hrms = hrms + cont(j)*(v - hmean)**2
      hskew = hskew + cont(j)*(v - hmean)**3
      hkur = hkur + cont(j)*(v - hmean)**4
      enddo
      hrms = sqrt(hrms/fn)
      hskew = hskew/fn/hrms**3
      hkur = hkur/fn/hrms**4
      write(*,*)' '
      write(*,*)' HITS 2, mean rms  skewness kurtosis',
     1 ii,hmean,hrms,hskew,hkur
      write(*,*)' '
      vhskew1(1) = hskew
      vhkur1(1) = hkur

c234567
      call hunpak(12000,cont,' ',0) ! projection folded
      frms = 0
      fmean = 0.
      do j = 1,230
      fj = (j-1)*4  
      temp = cont(j+1) - cont(j) ! first diff uncorrected folded distr.
      frms = frms + temp**2
      fmean = fmean  + temp
      call hfill(12100,fj,0.,temp)
      call hfill(12200,temp,0.,1.)
      if(j.lt.229)then
      temp = cont(j) + cont(j+2) -2.*cont(j+1)
      call hfill(12300,fj,0.,temp)
      endif 
      enddo
      fmean = fmean/230.
      frms = sqrt(frms/230.)

      write(*,*)' '
      write(*,*)' First Diff.  mean ,rms ',fmean,frms
      write(*,*)'  '
      vfrms(1) = frms


c--------------------------------------------
*     call hreset(100,' ')
*     call hreset(101,' ')
*     call hreset(102,' ')
*     call hreset(103,' ')
*     call hreset(104,' ')

*********************************      enddo  ! loop over files
      go to 111
      call hprint(200)
      call hprint(201)
      call hprint(202)
      call hprint(203)
      call hprint(204)
      call hprint(205)
      call hprint(206)
      call hprint(207)

      call hprint(301)
      call hprint(302)
      call hprint(303)
      call hprint(304)
      call hprint(305)
      call hprint(306)
      call hprint(307)
 111  continue

*     periodogram / correlogram
      do ii =1,12
*     call hprint(9200+ii)
      enddo


c      call hrput(0,'spillstat3.dat','n')
      stop
      end
quit

call spillstat4

quit


*h/file 1 spillstat3.dat
*fort/file 20 spill4.ps
* spill4.kumac
* reads output from spillstat4.f
* was spillstat4.kumac
meta 20 -111
 igset mtype 20
 igset mscf 1.5
set yhti 100.
set ygti 1.0
set xmgl 3.
set xlab 1.00
opt stat
*set ndvy 50.
igset lwid 2.
set hwid 5
opt date
set xwin 2.5
* helvetica
igset txfp -40 
set lfon -40
set vfon -40
set gfon -40
h/list

go to aaa
zone 2 2
title ' Spill statistics (from 1.5 - 5.5 secs spill time, 1 ms interval)  '
h/pl 200
atitle ' mean (triggers/ms) ' 'events'
h/pl 201
atitle ' RMS/sqrt(mean) ' 
h/pl 202
atitle ' skewness (0 for Gaussian) ' 'events'
h/pl 203
atitle ' kurtosis (3 for Gaussian)' ' '

title ' noise fractions '
zone 2 2
h/pl 204
atitle ' total noise below 40 Hz ' ' events '
h/pl 205
atitle ' 100 Hz noise fraction '
h/pl 206
atitle ' total noise fraction '
h/pl 207
atitle ' High spill % '

zone 2 2
opt nstat
set xlab 1.50

h/pl 301

atitle ' RMS/sqrt(mean) ' 'skewness'
h/pl 302
atitle ' RMS/sqrt(mean) ' 'kurtosis'
h/pl 303
atitle ' skewness ' 'kurtosis'
zone 2 2
h/pl 304
atitle ' RMS/sqrt(mean) ' ' low noise fraction'
h/pl 305
atitle ' RMS/sqrt(mean) ' ' 100 Hz fraction'
h/pl 306
atitle ' RMS/sqrt(mean) ' ' total noise  fraction'
h/pl 307
atitle ' High folded spill % '

aaa:

set xlab 1.50
opt nstat
opt nfit
set hwid 3
opt grid

tempb = ibur(1)
tempr = irun(1)
temp6 = ispi(1)
tempn = inoi(1)
tempf = ispf(1)
t100  = n100(1)
tfreq = ifreq(1)
trms  = irms(1)
trms1 = irms1(1) 
temp = scal(1)
tempx = imx5(1)
temp5 = imn5(1)


 igset mtype 20
 igset mscf 0.1
 opt nstat
 opt nfit 
 temp = scal(1)
zone 1 1
title ' corrected spill (2D spill cut 10, hit cut 0)'
id = 10000 
h/pl [id]
itx 100 5900 '  ' [tempr] '  ' [tempb] 
itx 400 5900 'TIME ' [temp]
temp = vcrem(1)
itx 650 5900 '% removed ' [temp]
atitle 'Fig. 1           trigger time modulo(923.99xxxx) (25 ns) ' ' trigger time (ms) '

zone 1 4
title ' Spill - full range - 5 ms bins '
id = 110 
h/pl [id]
itx 2000 800 '  ' [tempr] '  ' [tempb] 
itx 3000 800 ' ' [temp]
set txci 1
if [tempx] >  800 then
set txci 2
endif
itx 4000 800 'MAX ' [tempx]
set txci 1
if [temp5] <  100 then
set txci 2
endif
itx 4750 800 'MIN ' [temp5]
set txci 1
*itx 6100  900 'excess '  
*itx 6300 650 [temp6] ' %'
itx 6100 400 'rms 1 ms '
itx 6300 150 [trms1] 
set plci 7
set lwid 5
dline 1000 6000 600 600
set plci 1
set lwid 1
atitle 'Fig. 2      spill (ms) ' ' Triggers / 5 ms'
title '  corrected folded spill '
*opt grid
*id =  11000 + [i]
*h/pl [id]
*atitle ' corrected folded spill '

set hwid 3
title ' folded spill '
id =  12000 
opt grid
h/pl [id]
set plci 7
set lwid 5
*dline 0 924 500 500
set plci 1
set lwid 1
*itx 940 500 'excess ' 
*itx 960 250 [tempf] ' %'
if [trms] >  100 then
set txci 2
endif
itx 940 00  'rms ' [trms] 
set txci 1
atitle 'Fig.3      Uncorrected folded spill (1.2 - 2.2 secs) ' ' Triggers'

opt grid
set hwid 3
*zone 1 2
title ' Spill Frequency Analysis from Periodogram/Autocorrelation  (1.5 - 2.5 secs) '
id =  8000 
h/pl [id]
if [tempn] >  45 then
set txci 2
endif
*itx 200 2000 [tempn] ' %'
 itx 260 2500  'noise '
 itx 260  000 [tempn] ' %'
set txci 1
atitle 'Fig. 4      Frequency (Hz)' ' Amplitude'

set hwid 2
title ' Autocorrelation'
 
id =  9200 
h/pl [id](001.:450.)
if [t100] >  8000 then
set txci 2
endif
 itx 455  0.20  'Max amp '
 itx 455  0 [t100] 
 itx 455 -0.2 [tfreq] ' Hz' 
set txci 1
atitle 'Fig. 5              Time (ms)' ' Autocorrelation' 
opt ngrid
set hwid 5

igset mtype 1
igset mscf .5 

opt stat
opt fit 

go to zzz
title ' Spill:  counts / 5 ms (1.5 - 5.5 secs).       First difference for 1 ms '
zone 2 2
set csiz 0.35
id = 9000 
h/fit [id] G
set plci 7
set lwid 5
dline 600. 600. 0. 40.
set plci 1
set lwid 1
atitle '  events/5 ms ' ' entries'


title ' First difference of spill / ms  (1.5 - 5.5 secs)'
id = 9100 
h/fit [id] G
atitle '  events / 1 ms  difference ' ' entries'


 igset mtype 20
 igset mscf 0.1
 opt nstat
 opt nfit 
 temp = scal(1)
zone 1 1
title ' corrected spill '
id = 10000 
h/pl [id]
itx 350 5900 '  ' [tempr] '  ' [tempb] 
itx 650 5900 'TIME ' [temp]
atitle 'Fig. 5           trigger time modulo(923.99xxxx) (25 ns) ' ' trigger time (ms) '


zzz:



opt liny
opt nstat
opt nfit
set fwid 1
set plci 1 
opt grid
opt tic
title  '  Hits GTK(1) vs Spill Time, Folded Spill Time '
zone 1 2
 igset mtype 20
 igset mscf 0.11
 igset mscf 0.10
h/pl 20000
dline 2000 2000 0 250
dline 5000 5000 0 250
atitle 'Fig. 6         Spill Time (ms) ' ' Hits '
h/pl 21000
atitle 'Fig. 7   Folded  Spill Time (25 ns) ' ' Hits '




 igset mtype 20
 igset mscf 0.1
 opt nstat
 opt nfit 
temp = vhskew(1)
temp1 = vhkur(1)
opt stat
set plci 2
set fwid 5
set csiz 0.40
title ' Hits GTK(1) '
zone 2 2
h/fit  7000  G
if [temp] > 1 then
set txci 2
endif
itx 100 10000 ' Skewness ' [temp]
itx 100 7500 ' Kurtosis ' [temp1]
set txci 1
atitle 'Fig. 8a           hits ' 'entries'
opt logy
h/pl 7000
atitle 'Fig. 8b           hits ' 'entries'
opt liny 
temp = vhskew1(1)
temp1 = vhkur1(1)
h/fit 7100 G
itx 100 8000 ' Skewness ' [temp]
itx 100 6000 ' Kurtosis ' [temp1]
atitle 'Fig 8c   hits (2.0 -5.0 secs)'  'entries'
opt logy
h/pl 7100  
atitle 'Fig 8d   hits (2.0 -5.0 secs)'   'entries'
*--------------------
go to xyz
opt liny
opt nstat
opt nfit
set fwid 1
set plci 1 
opt grid
opt tic
title  '  Folded Spill Time  vs Hits GTK(1)'
zone 1 2
 igset mtype 20
 igset mscf 0.11
 igset mscf 0.10
h/pl 20000
dline 2000 2000 0 250
dline 5000 5000 0 250
atitle 'Fig. 7         Spill Time (ms) ' ' Hits '
h/pl 21000
atitle 'Fig. 8   Folded  Spill Time (25 ns) ' ' Hits '
xyz:
*----------------------------------------
opt ngrid
opt ntic
opt stat
opt fit
set csiz 0.3

go to www
title ' fist diff. uncorr.folded time '
h/pl 12100
atitle 'Fig. 9         uncorr. folded time ' ' first diff '
h/fit 12200 G
temp = vfrms(1)
if [temp] > 100 then
set txci 2
endif
itx 400 5 'RMS  ' [temp]
set txci 1
atitle 'Fig. 10         first diff folded spill time '  ' events '
opt nstat
opt nfit
*title ' Second difference '
*h/pl 12300
*atitle '  Folded spill time (25 ns) ' ' Second difference '
www:

opt nstat
opt nfit
zone 1 3
title 'Comparison of GTK(1) hits /5ms  with triggers /5 ms  '
opt grid
opt logy
h/pl 31110  
atitle 'Fig 11a             spill time (ms)'   ' hits  triggers'
set dmod 1
h/pl 30110 s
*atitle ' spill time (ms)'  ' triggers'
opt liny
h/pl 32110   
atitle 'Fig. 11b             spill time (ms)'   '  hits / triggers'
h/pl 33110   
atitle 'Fig. 11c             spill time (ms)'   '  triggers / hits'
opt ngrid

enddo


quit