macro dfact field=perp ttype=h fit=data
read field
fittype=data
fittype=data
pmom=4.723
* Lower and upper region of W to match background and data
fwlo=0.6
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]
btot_charge=0
btot_eff=0
first=4
last=[nrtot]
diff=[last]-[first]+1
do nr=[first],[last]
nrun=vrun([nr])
if ($FEXIST(t.out)) then
sh rm t.out
endif
sh tclsh ../get_info_run.tcl [nrun] > t.out
sh cat t.out
vec/del ch,t,f
vec/read ch,t,f t.out
vec/write ch,t,f
bcharge=ch(1)
message charge = [bcharge] mC
btime=t(1)
message time = [btime] sec
bcurrent=[bcharge]/[btime]*1000.
message ave current = [bcurrent] nA
bfideff=f(1)
bfideff=[bfideff]/100.
message Tracking eff = [bfideff]
btot_charge=[btot_charge]+[bcharge]
btot_eff=[btot_eff]+[bfideff]
sh rm  t.out
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
bavefideff=[btot_eff]/[diff]
bcf=1./[btot_charge]/[bavefideff]
message            Total charge = [btot_charge] 
message            average fideff= [bavefideff]
*
* calc dilution factor 
*
ntheta=2
if ( [field] .eq. 'para') then
      vec/create vthlo([ntheta]) r 11.5    13.0875
      vec/create vthhi([ntheta]) r 13.0875 15.0
endif
if ( [field] .eq. 'perp') then
      vec/create vthlo([ntheta]) r  11.5    12.7625
      vec/create vthhi([ntheta]) r  12.7625 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
     if ( [fittype] .eq. 'yield') then
	vec/read bgw,bgc,bgerr [field]_yield_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
     if ( [fittype] .eq. 'yield') then
	vec/read bgw,bgc,bgerr [field]_yield_bg_43396_th[nth].dat
     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])
if ([field] .eq. 'perp') then
ntu/cut $11 abs(hsyptar)<.04.and.abs(hsxptar+.05)<.12
endif
if ([field] .eq. 'para') then
ntu/cut $11 abs(hsyptar)<.04.and.abs(hsxptar)<.12
endif
ntu/cut $2 $11.and.[thlo]<hstheta<=[thhi]
nbin=50
wlo=.6
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.
vec/create errcon(1) r 
hist/fit 1011([fwlo]:[fwhi]) fit_back_data.f SQ 1  fcon ! ! ! errcon
endif
if ([fittype] .eq. 'yield') then
temp=.8/[bcf]
vec/create fcon(1) r [temp]
vec/create errcon(1) r 
vec/create evw([nbin]) r [nbin]*.001
*hplot/errors vw tot evw etot 200 24 .1 W
hist/fit 1011([fwlo]:[fwhi]) fit_back_data.f SQ 1  fcon ! ! ! errcon
endif

vec/write fcon,errcon
con=fcon(1) 
*
if ([fittype] .eq. 'data') then
con=fcon(1)
sigma yy=bgc*[con]
endif
if ([fittype] .eq. 'yield') then
con=fcon(1)
temp=[con]/.8*[bcf]
etemp=errcon(1)/.8*[bcf]
message Fit constant is [temp] +- [etemp] relative to 0.8*[charge ratio]
sigma yy=bgc*[con]
sigma yyerr=bgerr*[con]
set plci 2.
*hplot/errors vw yy evw yyerr 200 24 .1 
set plci 1.
vec/write  bgw,yy,yyerr,tot,etot ' ' 5(f15.5,1x)
*wait
endif
SET PLCI 4 
pline [nbin] 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
if ( [fittype] .eq. 'yield') then
sigma edf=sqrt(yy/tot/tot*etot*yy/tot/tot*etot+yyerr*yyerr/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