macro asym_plot wtype=elas field=perp ttype=h
opt fit
set ltyp 1
set hcol 1.
set dmod 1
vec/del *
title_gl '  '
pmom=4.723
start_theta=11.0
if ([field] .eq. 'perp') then
start_theta=11.5
endif
message field = [field]
message target = [ttype]
message W region = [wtype]
*sh tclsh get_pol.tcl [field][ttype] [pmom]
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].dat
  o1name=top_th1_[field][ttype]_[wtype].dat
  o11name=top_th11_[field][ttype]_[wtype].dat
  o2name=top_th2_[field][ttype]_[wtype].dat
  o12name=top_th12_[field][ttype]_[wtype].dat
  o3name=top_th3_[field][ttype]_[wtype].dat
  o4name=top_th4_[field][ttype]_[wtype].dat
  w1name=top_w1_[field][ttype]_[wtype].dat
  w2name=top_w2_[field][ttype]_[wtype].dat
  w3name=top_w3_[field][ttype]_[wtype].dat
  w4name=top_w4_[field][ttype]_[wtype].dat
endif
if ( [ntar] .eq. 2) then
  fname=bot_[field][ttype]_[pmom].out
  oname=bot_[field][ttype]_[wtype].dat
  o1name=bot_th1_[field][ttype]_[wtype].dat
  o11name=bot_th11_[field][ttype]_[wtype].dat
  o2name=bot_th2_[field][ttype]_[wtype].dat
  o12name=bot_th12_[field][ttype]_[wtype].dat
  o3name=bot_th3_[field][ttype]_[wtype].dat
  o4name=bot_th4_[field][ttype]_[wtype].dat
  w1name=bot_w1_[field][ttype]_[wtype].dat
  w2name=bot_w2_[field][ttype]_[wtype].dat
  w3name=bot_w3_[field][ttype]_[wtype].dat
  w4name=bot_w4_[field][ttype]_[wtype].dat
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
vec/create vwpeak([nrtot]) r
vec/create velas([nrtot]) r
vec/create vwsig([nrtot]) r
vec/create vdf([nrtot]) r
vec/create vth([nrtot]) r
vec/create vth1([nrtot]) r
vec/create vth11([nrtot]) r
vec/create vth2([nrtot]) r
vec/create vth12([nrtot]) r
vec/create vth3([nrtot]) r
vec/create vth4([nrtot]) r
vec/create vw1([nrtot]) r
vec/create vwd1([nrtot]) r
vec/create vwd2([nrtot]) r
vec/create vwd3([nrtot]) r
vec/create vwd4([nrtot]) r
vec/create vw2([nrtot]) r
vec/create vw3([nrtot]) r
vec/create vw4([nrtot]) r
vec/create vasy([nrtot]) r
vec/create vasyerr([nrtot]) r
vec/create vasyc([nrtot]) r
vec/create vasycerr([nrtot]) r
vec/create vasy1c([nrtot]) r
vec/create vasy1cerr([nrtot]) r
vec/create vasy11c([nrtot]) r
vec/create vasy11cerr([nrtot]) r
vec/create vasy2c([nrtot]) r
vec/create vasy2cerr([nrtot]) r
vec/create vasy12c([nrtot]) r
vec/create vasy12cerr([nrtot]) r
vec/create vasy3c([nrtot]) r
vec/create vasy3cerr([nrtot]) r
vec/create vasy4c([nrtot]) r
vec/create vasy4cerr([nrtot]) r
vec/create wasy1c([nrtot]) r
vec/create wasy1cerr([nrtot]) r
vec/create wasy2c([nrtot]) r
vec/create wasy2cerr([nrtot]) r
vec/create wasy3c([nrtot]) r
vec/create wasy3cerr([nrtot]) r
vec/create wasy4c([nrtot]) r
vec/create wasy4cerr([nrtot]) r
*
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])
*
nseg=1
fname=../ntup/hms[nrun].[nseg].hbook
chain -h[field]
while ($FEXIST([fname])) do
message add [fname] to chain
chain h[field] [fname]
nseg=[nseg]+1
fname=../ntup/hms[nrun].[nseg].hbook
endwhile
zone 1 1
*
* constants from run 43775 para 12C+he 
 if ([field] .eq. 'para') then
      vec/create fparf(4) r 0.13247E+00 -0.47851E+00 0.54514E+00 -0.18406E+00
      vec/create fpar1(4) r 0.34406E-01   -0.12129E+00    0.13397E+00   -0.43785E-01
      vec/create fpar2(4) r 0.52919E-01   -0.18597E+00    0.20898E+00   -0.72358E-01
      vec/create fpar3(4) r 0.25090E-01   -0.94794E-01    0.11137E+00   -0.37743E-01
      vec/create fpar4(4) r 0.15389E-01   -0.61279E-01    0.74606E-01   -0.24548E-01
      vec/create fpar11(4) r 0.88035E-01   -0.30962E+00    0.34558E+00   -0.11710E+00
      vec/create fpar12(4) r 0.40000E-01   -0.15455E+00    0.18441E+00   -0.61758E-01
      vec/create vthlo(4) r 11.25 12.25 12.8 13.5
      vec/create vthhi(4) r 12.25 12.8  13.5 15.0
      vec/create vth1lo(2) r 11.25 12.8 
      vec/create vth1hi(2) r 12.8  15.0
endif
* constants from run 43229 perp 12C+he
 if ([field] .eq. 'perp') then
      vec/create fparf(4) r 0.10957E+00 -0.39757E+00 0.45319E+00 -0.15035E+00
      vec/create fpar1(4) r  0.42405E-01   -0.14945E+00    0.16689E+00   -0.56459E-01
      vec/create fpar2(4) r  0.36657E-01   -0.13155E+00    0.14956E+00   -0.51008E-01
      vec/create fpar3(4) r  0.20426E-01   -0.78576E-01    0.93997E-01   -0.32178E-01
      vec/create fpar4(4) r  0.96519E-02   -0.37447E-01    0.43142E-01   -0.11260E-01
      vec/create fpar11(4) r 0.79430E-01   -0.28224E+00    0.31786E+00   -0.10799E+00
      vec/create fpar12(4) r 0.29975E-01   -0.11565E+00    0.13672E+00   -0.43289E-01
      vec/create vthlo(4) r 11.5 12.5  13.1 13.8
      vec/create vthhi(4) r  12.5 13.1  13.8 16.0
      vec/create vth1lo(2) r 11.5 13.1
      vec/create vth1hi(2) r  13.1 16.0
 endif
* calc dilution factor ( no theta cut)
*
nbin=200
wlo=.7
whi=1.100
step=([whi]-[wlo])/[nbin]
ttit= W  nrun=[nrun]
cdir //pawc
hist/del 0
cdir //h[field]
hist/del 1011
1dhisto 1011 [ttit] 80  .8 .88 
ntu/proj 1011 //h[field]/9010.w hsp>.01.and.0.8<hsshtrk/hsp<1.4
hrin 1011
sum1=$hinfo(1011,'sum')
1dhisto 1011 [ttit] [nbin]  [wlo] [whi] 
ntu/proj 1011 //h[field]/9010.w hsp>.01.and.0.8<hsshtrk/hsp<1.4
hist/plot 1011
vec/create fpar(4) r
vec/copy fparf fpar
vec/create vcon(1) r 1.
hist/fit 1011(.7:.85) fit_back.f SQ 1  vcon
vec/write vcon
con=vcon(1) 
*wait
****
*** This sets the W region to determine the df and asym

wflo=.9
wfhi=1.0
1dhisto 1012 [ttit] [nbin]  [wflo] [wfhi] 
ntu/proj 1012 //h[field]/9010.w hsp>.01.and.0.8<hsshtrk/hsp<1.4
hrin 1012
sum91=$hinfo(1012,'sum')
p0=fpar(1)
p1=fpar(2)
p2=fpar(3)
p3=fpar(4)
message sum91 = [sum91]
sigma x=array([nbin],[wlo]#[whi])
sigma y=[con]*([p0]+[p1]*x+[p2]*x*x+[p3]*x*x*x)
sigma sumy=sumv([con]*([p0]+[p1]*x+[p2]*x*x+[p3]*x*x*x))
*vec/write x,y,sumy ' ' 3e15.5
* get background from W between [wflo] and[wfhi]
bin9=int(([wflo]-[wlo])*[nbin]/([whi]-[wlo]))
bin1=int(([wfhi]-[wlo])*[nbin]/([whi]-[wlo]))
cnt9=sumy([bin9])
cnt1=sumy([bin1])
tot_back=[cnt1]-[cnt9]
pline [nbin] x y
df=([sum91]-[tot_back])/[sum91]
elas=[sum91]-[tot_back]
message [bin9] [cnt9] counts at W = [wflo] 
message [bin1] [cnt1] counts at W = [wfhi]
message tot_back = [tot_back] tot=[sum91] elas=[elas]   d.f. = [df]
vec/input vdf([nr]) [df] 
vec/input velas([nr]) [elas] 
*
if ( [sum1] .ne. 0) then
vec/create x([nbin]) r
vec/create y([nbin]) r
hist/get_vect/contents 1011 y
hist/get_vect/abscissa 1011 x
sigma maxc=vmax(y)
sigma maxw=lvmax(y)
wc=maxc(1)
wpeak=[wlo]+[step]*maxw(1)
wwid=.02
flo=[wpeak]-[wwid]
fhi=[wpeak]+[wwid]
vec/create par(3) r [wc] [wpeak] [wwid]
message fit with [wc] [wpeak] [wwid] 
set fcol 2.
hist/fit 1011([flo]:[fhi]) g SQ 3  par
vec/write par ' ' 4f11.5
vec/input vwpeak([nr]) par(2) 
vec/input vwsig([nr]) par(3) 
metafile 44 -111
*wait
*
* calc asym for all theta
*
1dhisto 102 ' theta_scat (deg) ' 100 10 16
ntu/proj 102 //h[field]/9010.hstheta hsp>.01.and.0.8<hsshtrk/hsp<1.4.and.[wflo]<W<[wfhi]
hrin 102
mean=$hinfo(102,'mean')
vec/input vth([nr]) [mean] 
1dhisto 1011 ' W h+' [nbin]  [wflo] [wfhi] 
ntu/proj 1011 //h[field]/9010.w hsp>.01.and.0.8<hsshtrk/hsp<1.4.and.helicity.ge.32768
hrin 1011
pct=$hinfo(1011,'sum')
1dhisto 1011 ' W h-' [nbin]  [wflo] [wfhi] 
ntu/proj 1011 //h[field]/9010.w hsp>.01.and.0.8<hsshtrk/hsp<1.4.and.helicity.lt.32768
hrin 1011
mct=$hinfo(1011,'sum')
asy=([mct]-[pct])/([pct]+[mct])
asyc=([asy]-[chasy])/(1-[asy]*[chasy])/[bpol]/[tpol]/[df]
sigma err=2*sqrt([pct]*[mct])/([pct]+[mct])/sqrt([pct]+[mct])
asyerr=err(1)
asyc_err=abs([asyerr]*(1-[chasy]*[chasy])/(1-[chasy]*[asy])/(1-[chasy]*[asy])/[bpol]/[tpol]/[df])
vec/input vasy([nr]) [asy] 
vec/input vasyc([nr]) [asyc] 
vec/input vasyerr([nr]) [asyerr] 
vec/input vasycerr([nr]) [asyc_err] 
message theta = [mean] asy=[asy]+-[asyerr] asy_cor=[asyc]+-[asyc_err]
*
*
do i=1,4
wflo=.90+([i]-1)*.025
wfhi=[wflo]+.025
ntu/cut $1 [wflo]<w<=[wfhi]
1dhisto 1012 [ttit] [nbin]  .7 1.2 
ntu/proj 1012 //h[field]/9010.w $1.and.hsp>.01.and.0.8<hsshtrk/hsp<1.4
hrin 1012
sum91=$hinfo(1012,'sum')
mean=$hinfo(1012,'mean')
message sum91 = [sum91]
p0=fpar(1)
p1=fpar(2)
p2=fpar(3)
p3=fpar(4)
message [i] [p0] [p1] [p2] [p3]
sigma x=array([nbin],[wlo]#[whi])
sigma y=[con]*([p0]+[p1]*x+[p2]*x*x+[p3]*x*x*x)
sigma sumy=sumv([con]*([p0]+[p1]*x+[p2]*x*x+[p3]*x*x*x))
*vec/write x,y,sumy ' ' 3e15.5
* get background from W between [wflo] and[wfhi]
bin9=int(([wflo]-[wlo])*[nbin]/([whi]-[wlo]))
bin1=int(([wfhi]-[wlo])*[nbin]/([whi]-[wlo]))
cnt9=sumy([bin9])
cnt1=sumy([bin1])
tot_back=[cnt1]-[cnt9]
pline [nbin] x y
df=([sum91]-[tot_back])/[sum91]
elas=[sum91]-[tot_back]
message [bin9] [cnt9] counts at W = [wflo] 
message [bin1] [cnt1] counts at W = [wfhi]
message tot_back = [tot_back] tot=[sum91] elas=[elas]   d.f. = [df]
1dhisto 1011 ' W h+' [nbin]  .7 1.2 
ntu/proj 1011 //h[field]/9010.w $1.and.hsp>.01.and.0.8<hsshtrk/hsp<1.4.and.helicity.ge.32768
hrin 1011
pct=$hinfo(1011,'sum')
1dhisto 1011 ' W h-' [nbin]  .7 1.2 
ntu/proj 1011 //h[field]/9010.w $1.and.hsp>.01.and.0.8<hsshtrk/hsp<1.4.and.helicity.lt.32768
hrin 1011
mct=$hinfo(1011,'sum')
asy=([mct]-[pct])/([pct]+[mct])
message asy  = [asy] [mct] [pct] 
asyc=([asy]-[chasy])/(1-[asy]*[chasy])/[bpol]/[tpol]/[df]
sigma err=2*sqrt([pct]*[mct])/([pct]+[mct])/sqrt([pct]+[mct])
asyerr=err(1)
asyc_err=abs([asyerr]*(1-[chasy]*[chasy])/(1-[chasy]*[asy])/(1-[chasy]*[asy])/[bpol]/[tpol]/[df])
vec/input vw[i]([nr]) [mean] 
vec/input vwd[i]([nr]) [df] 
vec/input wasy[i]c([nr]) [asyc] 
vec/input wasy[i]cerr([nr]) [asyc_err] 
message W = [mean] asy=[asy]+-[asyerr] asy_cor=[asyc]+-[asyc_err] df =  [df]
enddo
*
* calc asym for theta bins
*
zone 2 2
do nth=1,4
thlo=vthlo([nth])
thhi=vthhi([nth])
ntu/cut $1 [thlo]<hstheta<=[thhi]
1dhisto 1011 [ttit] [nbin]  [wlo] [whi] 
ntu/proj 1011 //h[field]/9010.w $1.and.hsp>.01.and.0.8<hsshtrk/hsp<1.4
hist/plot 1011
vec/create fpar(4) r
vec/copy fpar[nth] fpar
vec/write fpar ' ' 4(1x,e15.5)
vec/create vcon(1) r 1.
hist/fit 1011(.7:.85) fit_back.f SQ 1  vcon
vec/write vcon
con=vcon(1) 
wflo=.90
wfhi=1.0
1dhisto 1012 [ttit] [nbin]  [wflo] [wfhi] 
ntu/proj 1012 //h[field]/9010.w $1.and.hsp>.01.and.0.8<hsshtrk/hsp<1.4
hrin 1012
sum91=$hinfo(1012,'sum')
message sum91 = [sum91]
p0=fpar(1)
p1=fpar(2)
p2=fpar(3)
p3=fpar(4)
message nth = [nth] [p0] [p1] [p2] [p3]
sigma x=array([nbin],[wlo]#[whi])
sigma y=[con]*([p0]+[p1]*x+[p2]*x*x+[p3]*x*x*x)
sigma sumy=sumv([con]*([p0]+[p1]*x+[p2]*x*x+[p3]*x*x*x))
*vec/write x,y,sumy ' ' 3e15.5
* get background from W between [wflo] and[wfhi]
bin9=int(([wflo]-[wlo])*[nbin]/([whi]-[wlo]))
bin1=int(([wfhi]-[wlo])*[nbin]/([whi]-[wlo]))
cnt9=sumy([bin9])
cnt1=sumy([bin1])
tot_back=[cnt1]-[cnt9]
pline [nbin] x y
df=([sum91]-[tot_back])/[sum91]
1dhisto 102 ' theta_scat (deg) ' 100 10 16
ntu/proj 102 //h[field]/9010.hstheta $1.and.hsp>.01.and.hsp>.01.and.0.8<hsshtrk/hsp<1.4.and.[wflo]<W<[wfhi]
hrin 102
mean=$hinfo(102,'mean')
vec/input vth[nth]([nr]) [mean] 
1dhisto 1011 ' W h+' [nbin]  [wflo] [wfhi] 
ntu/proj 1011 //h[field]/9010.w $1.and.hsp>.01.and.hsp>.01.and.0.8<hsshtrk/hsp<1.4.and.helicity.ge.32768
hrin 1011
pct=$hinfo(1011,'sum')
1dhisto 1011 ' W h-' [nbin]  [wflo] [wfhi] 
ntu/proj 1011 //h[field]/9010.w $1.and.hsp>.01.and.hsp>.01.and.0.8<hsshtrk/hsp<1.4.and.helicity.lt.32768
hrin 1011
mct=$hinfo(1011,'sum')
asy=([mct]-[pct])/([pct]+[mct])
sigma err=2*sqrt([pct]*[mct])/([pct]+[mct])/sqrt([pct]+[mct])
asyerr=err(1)
asyc=([asy]-[chasy])/(1-[asy]*[chasy])/[bpol]/[tpol]/[df]
asyc_err=abs([asyerr]*(1-[chasy]*[chasy])/(1-[chasy]*[asy])/(1-[chasy]*[asy])/[bpol]/[tpol]/[df])
vec/input vasy[nth]c([nr]) [asyc] 
vec/input vasy[nth]cerr([nr]) [asyc_err] 
message theta = [mean] asy=[asy]+-[asyerr] asy_cor=[asyc]+-[asyc_err] df = [df]
enddo
*
* calc asym for theta bins
*
zone 2 2
do nth=1,2
thlo=vth1lo([nth])
thhi=vth1hi([nth])
ntu/cut $1 [thlo]<hstheta<=[thhi]
1dhisto 1011 [ttit] [nbin]  [wlo] [whi] 
ntu/proj 1011 //h[field]/9010.w $1.and.hsp>.01.and.0.8<hsshtrk/hsp<1.4
hist/plot 1011
vec/create fpar(4) r
vec/copy fpar1[nth] fpar
vec/write fpar ' ' 4(1x,e15.5)
vec/create vcon(1) r 1.
hist/fit 1011(.7:.85) fit_back.f SQ 1  vcon
vec/write vcon
con=vcon(1) 
wflo=.90
wfhi=1.0
1dhisto 1012 [ttit] [nbin]  [wflo] [wfhi] 
ntu/proj 1012 //h[field]/9010.w $1.and.hsp>.01.and.0.8<hsshtrk/hsp<1.4
hrin 1012
sum91=$hinfo(1012,'sum')
message sum91 = [sum91]
p0=fpar(1)
p1=fpar(2)
p2=fpar(3)
p3=fpar(4)
message nth = [nth] [p0] [p1] [p2] [p3]
sigma x=array([nbin],[wlo]#[whi])
sigma y=[con]*([p0]+[p1]*x+[p2]*x*x+[p3]*x*x*x)
sigma sumy=sumv([con]*([p0]+[p1]*x+[p2]*x*x+[p3]*x*x*x))
*vec/write x,y,sumy ' ' 3e15.5
* get background from W between [wflo] and[wfhi]
bin9=int(([wflo]-[wlo])*[nbin]/([whi]-[wlo]))
bin1=int(([wfhi]-[wlo])*[nbin]/([whi]-[wlo]))
cnt9=sumy([bin9])
cnt1=sumy([bin1])
tot_back=[cnt1]-[cnt9]
pline [nbin] x y
df=([sum91]-[tot_back])/[sum91]
1dhisto 102 ' theta_scat (deg) ' 100 10 16
ntu/proj 102 //h[field]/9010.hstheta $1.and.hsp>.01.and.hsp>.01.and.0.8<hsshtrk/hsp<1.4.and.[wflo]<W<[wfhi]
hrin 102
mean=$hinfo(102,'mean')
vec/input vth1[nth]([nr]) [mean] 
1dhisto 1011 ' W h+' [nbin]  [wflo] [wfhi] 
ntu/proj 1011 //h[field]/9010.w $1.and.hsp>.01.and.hsp>.01.and.0.8<hsshtrk/hsp<1.4.and.helicity.ge.32768
hrin 1011
pct=$hinfo(1011,'sum')
1dhisto 1011 ' W h-' [nbin]  [wflo] [wfhi] 
ntu/proj 1011 //h[field]/9010.w $1.and.hsp>.01.and.hsp>.01.and.0.8<hsshtrk/hsp<1.4.and.helicity.lt.32768
hrin 1011
mct=$hinfo(1011,'sum')
asy=([mct]-[pct])/([pct]+[mct])
sigma err=2*sqrt([pct]*[mct])/([pct]+[mct])/sqrt([pct]+[mct])
asyerr=err(1)
asyc=([asy]-[chasy])/(1-[asy]*[chasy])/[bpol]/[tpol]/[df]
asyc_err=abs([asyerr]*(1-[chasy]*[chasy])/(1-[chasy]*[asy])/(1-[chasy]*[asy])/[bpol]/[tpol]/[df])
vec/input vasy1[nth]c([nr]) [asyc] 
vec/input vasy1[nth]cerr([nr]) [asyc_err] 
message theta = [mean] asy=[asy]+-[asyerr] asy_cor=[asyc]+-[asyc_err] df = [df]
enddo
*
endif
chain -h[field]
*wait
enddo
***
vec/write vrun,vwpeak,vwsig,vtpol,vbpol,vdf,vth,vcasy,vasy,vasyerr,vasyc,vasycerr,charge,time,acur,velas ' '     (f7.1,2(',',f8.5),3(',',f5.3),',',f6.3,9(',',e12.5))
vec/write vrun,vwpeak,vwsig,vtpol,vbpol,vdf,vth,vcasy,vasy,vasyerr,vasyc,vasycerr,charge,time,acur,velas [oname] (f7.1,2(',',f8.5),3(',',f5.3),',',f6.3,9(',',e12.5))
vec/write vrun,vth1,vasy1c,vasy1cerr [o1name] (f7.1,',',f6.3,2(',',e12.5))
vec/write vrun,vth11,vasy11c,vasy11cerr [o11name] (f7.1,',',f6.3,2(',',e12.5))
vec/write vrun,vth2,vasy2c,vasy2cerr [o2name] (f7.1,',',f6.3,2(',',e12.5))
vec/write vrun,vth12,vasy12c,vasy12cerr [o12name] (f7.1,',',f6.3,2(',',e12.5))
vec/write vrun,vth3,vasy3c,vasy3cerr [o3name] (f7.1,',',f6.3,2(',',e12.5))
vec/write vrun,vth4,vasy4c,vasy4cerr [o4name] (f7.1,',',f6.3,2(',',e12.5))
vec/write vrun,vw1,wasy1c,wasy1cerr,vwd1 [w1name] (f7.1,',',f6.3,3(',',e12.5))
vec/write vrun,vw2,wasy2c,wasy2cerr,vwd2 [w2name] (f7.1,',',f6.3,3(',',e12.5))
vec/write vrun,vw3,wasy3c,wasy3cerr,vwd3 [w3name] (f7.1,',',f6.3,3(',',e12.5))
vec/write vrun,vw4,wasy4c,wasy4cerr,vwd4 [w4name] (f7.1,',',f6.3,3(',',e12.5))
enddo
close 44