macro dfact field=perp ttype=h fit=data
read field
fittype=data
pmom=4.723
* Lower and upper region of W to match background and data
fwlo=0.7
fwhi=.85
wtype=w_[fwlo]_[fwhi]
message field = [field]
message target = [ttype]
message W region = [wtype]
do ntar=1,2
if ( [ntar] .eq. 1) then
  name=top
endif
if ( [ntar] .eq. 2) then
  name=bot
endif
  fname=../[name]_[field][ttype]_[pmom].out
vec/del  vrun,vtar,vtpol,vbpol,vcasy,ptime,ntime,pcharge,ncharge 
vec/read vrun,vtar,vtpol,vbpol,vcasy,ptime,ntime,pcharge,ncharge [fname]
nrtot=$vlen(vrun,1)
message Looping for [nrtot] runs in file [fname]
*
*
filecase keep
chain -all
do nr=1,[nrtot]
*do nr=1,3
nrun=vrun([nr])
nseg=1
dir=/work/hallc/e93026/rss/pass3
fname=[dir]/ntup/hms[nrun].[nseg].hbook
message check [fname]
while ($FEXIST([fname])) do
message add [fname] to chain
chain all [fname]
nseg=[nseg]+1
fname=[dir]/ntup/hms[nrun].[nseg].hbook
endwhile
zone 1 2
enddo
*
* calc dilution factor 
*
ntheta=4
if ( [field] .eq. 'para') then
      vec/create vthlo([ntheta]) r 11.5   12.225  12.775 13.425
      vec/create vthhi([ntheta]) r 12.225 12.775  13.425 15.0
endif
if ( [field] .eq. 'perp') then
      vec/create vthlo([ntheta]) r  11.5    12.4875  13.0875 13.8125
      vec/create vthhi([ntheta]) r  12.4875 13.0875  13.8125 16.0
endif
do nth=1,[ntheta]
if ([field] .eq. 'para') then
   if ( [fit] .eq. 'data') then
     if ( [fittype] .eq. 'data') then
	vec/read bgw,bgc,bgerr [field]_bg_43799_th[nth].dat
	vec/write  bgw,bgc,bgerr ' ' 3(f15.5,1x)
     endif
   endif
endif
if ([field] .eq. 'perp') then
*  
   if ( [fit] .eq. 'data') then
     if ( [fittype] .eq. 'data') then
	vec/read bgw,bgc,bgerr [field]_bg_43396_th[nth].dat
	vec/write  bgw,bgc,bgerr ' ' 3(f15.5,1x)
     endif
   endif
endif
  oname=[name]_[field][ttype]_[wtype]_df_[fittype]_[fit]_th[nth].dat
  psname=[name]_[field][ttype]_[wtype]_df_[fittype]_[fit]_th[nth].ps
thlo=vthlo([nth])
thhi=vthhi([nth])
ntu/cut $2 [thlo]<hstheta<=[thhi]
nbin=40
wlo=.7
whi=1.100
step=([whi]-[wlo])/[nbin] 
cdir //pawc
hist/del 0
cdir //all
ttit= W [field] [name] [thlo] [thhi]
title_gl [ttit] 
1dhisto 1011 [ttit] [nbin]  [wlo] [whi] 
ntu/proj 1011 //all/9010.w $2.and.hsp>.01.and.0.8<hsshtrk/hsp.and.abs(hsdelta)<12.
hist/plot 1011
vec/create vw([nbin]) r
vec/create tot([nbin]) r
vec/create etot([nbin]) r
hist/get_vect/contents 1011 tot
hist/get_vect/error 1011 etot
hist/get_vect/abscissa 1011 vw
if ([fittype] .eq. 'data') then
vec/create fcon(1) r 1000.
hist/fit 1011([fwlo]:[fwhi]) fit_back_data.f SQ 1  fcon
endif

vec/write fcon
con=fcon(1) 
*
if ([fittype] .eq. 'data') then
con=fcon(1)
sigma yy=bgc*[con]
endif
SET PLCI 4 
pline 40 vw yy
SET PLCI 1
sigma df=1-yy/tot
if ( [fittype] .eq. 'data') then
sigma edf=sqrt(yy/tot/tot*etot*yy/tot/tot*etot+bgerr*bgerr*[con]*[con]/tot/tot)
endif
vec/create evw([nbin]) r [nbin]*.001
hplot/errors vw df evw edf 200 24 .1 W
vec/write vw,df,edf [oname] 3(f15.5,1x)
pict/print [psname] 
enddo
*
chain -all
enddo