field=para
read field
vec/del *
filecase keep
mcdir=df_qfs_scalein1.0_dec22
* Bottom perp
*mcname=hi[field]_43364_df_nh3_43364_60.2_43364_60.2_hi
* Top perp
mcname=hi[field]_43369_df_nh3_43369_57.5_43369_57.5_hi
*mcdir=../../pol_hms_single/outfiles/hbook
*mcname=[field]hi_nh3_pf55_large
mcfile=[mcdir]/[mcname].1.hbook
message Opening [mcfile]
nseg=1
message [mcfile]
chain -mc
while ($FEXIST([mcfile])) do
message add [mcfile] to chain
chain mc [mcfile]
nseg=[nseg]+1
mcfile=[mcdir]/[mcname].[nseg].hbook
endwhile
*
mcdir=../../pol_hms_single/outfiles/hbook
mcname=[field]hi_che_large
mcfile=[mcdir]/[mcname].1.hbook
nseg=1
message [mcfile]
chain -mc2
while ($FEXIST([mcfile])) do
message add [mcfile] to chain
chain mc2 [mcfile]
nseg=[nseg]+1
mcfile=[mcdir]/[mcname].[nseg].hbook
endwhile
*
zone 1 3
cut $1 numer_wt.and.a>1.and.abs(hsdelta)<8.
cut $2 numer_wt.and.abs(hsdelta)<8.
nbin=25
wlo=.6
whi=1.5
1dhist 1000 'W Nh3 background' [nbin] [wlo] [whi]
ntu/proj 1000 //mc/1.W $1
hist/plot 1000 
vec/create rate([nbin]) r
hist/get_vect/contents 1000 rate
1dhist 1000 'W CHe background' [nbin] [wlo] [whi]
ntu/proj 1000 //mc2/1.W $2
set hcol 2.
hist/plot 1000 s
vec/create crate([nbin]) r
hist/get_vect/contents 1000 crate
set hcol 1
1dhist 1000 'W Nh3 background' [nbin] [wlo] [whi]
ntu/proj 1000 //mc/1.W a>1.and.abs(hsdelta)<8.
h/plot 1000
vec/create wv([nbin]) r
vec/create wverr([nbin]) r [nbin]*0.0001
vec/create tot([nbin]) r
vec/create etot([nbin]) r
hist/get_vect/contents 1000 tot
hist/get_vect/error 1000 etot
hist/get_vect/abscissa 1000 wv
sigma ratio=rate/tot
sigma erate=ratio*sqrt(tot)
1dhist 1000 'W Che' [nbin] [wlo] [whi]
ntu/proj 1000 //mc2/1.W abs(hsdelta)<8.
set hcol 2.
h/plot 1000 s
set hcol 1.
vec/create btot([nbin]) r
vec/create ebtot([nbin]) r
hist/get_vect/contents 1000 btot
hist/get_vect/error 1000 ebtot
sigma cratio=crate/btot
sigma ecrate=cratio*sqrt(btot)
sigma rat=crate/rate
sigma errrat=rat*sqrt(ecrate*ecrate/crate/crate+erate*erate/rate/rate)
hplot/errors wv rate wverr erate [nbin] 24 .1 W
set pcol 2.
hplot/errors wv crate wverr ecrate [nbin] 20 .1 
set pcol 1.
wait
zone 1 1
hplot/errors wv rat wverr errrat [nbin] 24 .1 W
sigma wserr=vsum(1/errrat/errrat)
sigma werr=sqrt(1./wserr)
sigma wave=vsum(rat/errrat/errrat)
wt=wave(1)/wserr(1)
we=werr(1)
Message [field] average ratio = [wt] +- [we]
vec/create par(2)
vec/fit wv rat errrat p1 s 2 par 
vec/write par
vec/write wv,rat,errrat [field]_mc_ratio_15N_che.dat 3(f15.5,1x)