macro asym_plot wtype=elas field=perp 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]
*sh tclsh get_pol.tcl [field][ttype] [pmom]
* input the nitrogen asym from region 0.7 < W < .85
if ( [field] .eq. 'para') then
asymn = 0.003
endif
if ( [field] .eq. 'perp') then
asymn = 0.008
endif
* read in 15N dilution factor written by nitrogen_df.kumac
vec/read nvw,ndf [field]_15n_df.dat 
*
FOR/FILE 44 [field][ttype]_[wtype].ps
do ntar=1,2
if ( [ntar] .eq. 1) then
  fname=top_[field][ttype]_[pmom].out
  oname=top_[field][ttype]_[wtype]_count_asym_pass3.dat
  scal = 1.019
endif
if ( [ntar] .eq. 2) then
  fname=bot_[field][ttype]_[pmom].out
  oname=bot_[field][ttype]_[wtype]_count_asym_pass3.dat
  scal = 1.018
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
* 15N polarization
sigma vnpol=-(0.136*vtpol-.183*vtpol*vtpol+.335*vtpol*vtpol*vtpol)
*
message Looping for [nrtot] runs in file [fname]
*
*
do nr=1,[nrtot]
zone 1 1
*do nr=31,31
*
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
zone 1 2
*
nbin=40
if ([nr] .eq. 1) then
vec/create wsum([nbin]) r
vec/create esum([nbin]) r
vec/create ncsum([nbin]) r
endif
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
txt = Run [nrun]
1dhisto 1011 [txt] [nbin]  .7 1.1 
ntu/proj 1011 //h[field]/9010.w hsp>.01.and.0.8<hsshtrk/hsp<1.4.and.helicity.ge.32768.and.abs(hsdelta)<8.
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]  .7 1.1 
ntu/proj 1011 //h[field]/9010.w hsp>.01.and.0.8<hsshtrk/hsp<1.4.and.helicity.lt.32768.and.abs(hsdelta)<8.
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 ncorr=ndf*[npol]/[tpol]*[asymn]
sigma ncsum = ncsum + ncorr
sigma wsum = wsum + asyc/aerr/aerr
sigma esum = esum + 1./aerr/aerr
hplot/errors vw asyc vwerr aerr [nbin] 24 .3 W
xtxt = W ( GeV)  run [nrun]
ytxt = Count asymmetry
atitle [xtxt] [ytxt] ! 220 
*wait
else
  message [nrun] does not exist
endif
enddo
sigma werr=sqrt(1/esum)
sigma wave=wsum/esum
sigma ncsum=ncsum/[nrtot]
zone 1 2
hplot/errors vw wave vwerr werr [nbin] 24 .3 W
xtxt = W ( GeV)  [field] 
ytxt = Count asymmetry
atitle [xtxt] [ytxt] ! 220 
vec/write vw,wave,werr,ncsum [oname] 4(f10.6,1x)
*wait
*
close 44
enddo
***