macro asym_plot wtype=elas field=para ttype=h
read field
filecase keep
opt fit
set ltyp 1
set hcol 1.
set dmod 1
vec/del *
title_gl '  '
pmom=4.723
message field = [field]
message target = [ttype]
message W region = [wtype]
*
FOR/FILE 44 [field][ttype]_[wtype].ps
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 ntar=1,2
if ( [ntar] .eq. 1) then
  fname=../top_[field][ttype]_[pmom].out
  oname1=top_[field][ttype]_[wtype]_count_asym_pass3
  scal = 1.019
endif
if ( [ntar] .eq. 2) then
  fname=../bot_[field][ttype]_[pmom].out
  oname1=bot_[field][ttype]_[wtype]_count_asym_pass3
  scal = 1.018
  if ( [field] .eq. 'para') then
	scal = [scal]/1.18
  endif
endif
vec/del  vrun,vtar,vtpol,vbpol,vcasy,ptime,ntime,pcharge,ncharge 
vec/read vrun,vtar,vtpol,vbpol,vcasy,ptime,ntime,pcharge,ncharge [fname]
sigma charge=pcharge+ncharge
sigma time=ptime+ntime
sigma acur=charge/time*1000
nrtot=$vlen(vrun,1)
sigma vbpol=vbpol/100
*
message Looping for [nrtot] runs in file [fname]
*
*
zone 2 2
do nr=1,[nrtot]
*do nr=1,3
*
nrun=vrun([nr])
tpol=vtpol([nr])
bpol=vbpol([nr])
tar=vtar([nr])
chasy=vcasy([nr])
npol=vnpol([nr])
*
nseg=1
ddir=/work/hallc/e93026/rss/pass3
fname=[ddir]/ntup/hms[nrun].[nseg].hbook
chain -h[field]
if ( $FEXIST([fname]) ) then
while ($FEXIST([fname])) do
message add [fname] to chain
chain h[field] [fname]
nseg=[nseg]+1
fname=[ddir]/ntup/hms[nrun].[nseg].hbook
endwhile
*
nbin=50
if ([nr] .eq. 1) then
zone 2 2
do i=1,[ntheta]
vec/create wsum[i]([nbin]) r
vec/create esum[i]([nbin]) r
vec/create ncsum[i]([nbin]) r
enddo
endif
do nth=1,[ntheta]
vec/create vw([nbin]) r
vec/create vwerr([nbin]) r [nbin]*.01
vec/create pct([nbin]) r
vec/create perr([nbin]) r
vec/create mct([nbin]) r
vec/create merr([nbin]) r
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 $1 [thlo]<hstheta<=[thhi].and.$11
txt = Run [nrun]
1dhisto 1011 [txt] [nbin]  .6 1.1 
ntu/proj 1011 //h[field]/9010.w $1.and.hsp>.01.and.0.8<hsshtrk/hsp.and.helicity.ge.32768.and.abs(hsdelta)<12.
hrin 1011
hist/GET_VECT/abscissa 1011 vw 
hist/GET_VECT/CONTENTS 1011 pct
hist/GET_VECT/errors 1011 perr
1dhisto 1011 ' W h-' [nbin]  .6 1.1 
ntu/proj 1011 //h[field]/9010.w $1.and.hsp>.01.and.0.8<hsshtrk/hsp.and.helicity.lt.32768.and.abs(hsdelta)<12.
set hcol 2.
hrin 1011  
set hcol 1.
GET_VECT/CONTENTS 1011 mct
GET_VECT/errors 1011 merr
sigma asy=(mct-pct)/(pct+mct)
sigma asyc=-(asy-[chasy])/(1-asy*[chasy])/[bpol]/[tpol]/[scal]
sigma err=2*sqrt(pct*mct)/(pct+mct)/sqrt(pct+mct)
sigma aerr=abs(err*(1-[chasy]*[chasy])/(1-[chasy]*asy)/(1-[chasy]*asy)/[bpol]/[tpol]/[scal])
sigma wsum[nth] = wsum[nth] + asyc/aerr/aerr
sigma esum[nth] = esum[nth] + 1./aerr/aerr
hplot/errors vw asyc vwerr aerr [nbin] 24 .3 W
xtxt = W ( GeV)  run [nrun] [thlo] theta [thhi]
ytxt = Count asymmetry
atitle [xtxt] [ytxt] ! 220
enddo 
*wait
else
  message [nrun] does not exist
endif
enddo
zone 2 2
do nth=1,[ntheta]
sigma werr=sqrt(1/esum[nth])
sigma wave=wsum[nth]/esum[nth]
hplot/errors vw wave vwerr werr [nbin] 24 .3 W
xtxt = W ( GeV)  [field] theta [nth]
ytxt = Count asymmetry
atitle [xtxt] [ytxt] ! 220 
oname=[oname1]_th[nth].dat
vec/write vw,wave,werr [oname] 3(f10.6,1x)
enddo
*wait
*
enddo
***