macro dfact field=perp ttype=h fit=data
read field
fittype=g+p2
fittype=data
fittype=yield
pmom=4.723
* Lower and upper region of W to match background and data
fwlo=0.7
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
  oname=[name]_[field][ttype]_[wtype]_df_[fittype]_[fit].dat
  psname=[name]_[field][ttype]_[wtype]_df_[fittype]_[fit].ps
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
tot_charge=0
tot_eff=0
*do nr=1,[nrtot]
do nr=1,3
nrun=vrun([nr])
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
charge=ch(1)
time=t(1)
current=[charge]/[time]*1000.
fideff=f(1)/100.
sh rm t.out
message Total charge = [charge] mC
message Total time = [time] sec
message ave current = [current] nA
message Tracking eff = [fideff]
tot_charge=[tot_charge]+[charge]
tot_eff=[tot_eff]+[fideff]
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
nr=[nr]-1
ave_fideff=[tot_eff]/[nr]
cf=1./[tot_charge]/[ave_fideff]
message Total charge = [tot_charge] and averarge fideff= [ave_fideff]
*
if ([field] .eq. 'para') then
* constants from run 43775 para 12C+he 
      vec/create fparf(4) r 0.13247E+00 -0.47851E+00 0.54514E+00 -0.18406E+00
   if ( [fit] .eq. 'mc') then
* fti gaussian + p4 + p5*w + p6*w*w
      vec/create fit2(6) r 0.18126E+01    0.10843E+01    0.13257E+00    0.17292E-01    0.21076E-01    0.23168E-01
   endif
   if ( [fit] .eq. 'data') then
     if ( [fittype] .eq. 'g+p2') then
* fti gaussian + p4 + p5*w + p6*w*w
* fit runs  43775 43776 43799
      vec/create fit2(6) r  0.16197E+05    0.10906E+01    0.13425E+00   -0.54608E+03    0.24145E+03    0.12846E+04
* fit one run 43775
*      vec/create fit2(6) r 0.52741E+04    0.10942E+01    0.13876E+00    0.81637E+01    0.43012E+02    0.43515E+02
     endif
     if ( [fittype] .eq. 'data') then
	vec/read bgw,bgc,bgerr [field]_bg_43799.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.dat
	vec/write  bgw,bgc,bgerr ' ' 3(f15.5,1x)
        vec/create ypar(2) r 1.281 -0.044
     endif
   endif
endif
if ([field] .eq. 'perp') then
*  
      vec/create fparf(4) r 4*0.0
   if ( [fit] .eq. 'mc') then
     if ( [fittype] .eq. 'g+p2') then
* fti gaussian + p4 + p5*w + p6*w*w
      vec/create fit2(6) r 0.13274E+01    0.10774E+01    0.14874E+00    0.34963E-02    0.16148E-01    0.30040E-01
     endif
   endif
   if ( [fit] .eq. 'data') then
     if ( [fittype] .eq. 'g+p2') then
* fti gaussian + p4 + p5*w + p6*w*w
* fit runs  43227 43230 43285 43324 43395 43396
      vec/create fit2(6) r  0.26312E+05    0.11435E+01    0.18989E+00    0.62404E+04   -0.42061E+04   -0.88112E+04
* fit one run 43407
*      vec/create fit2(6) r 0.27566E+04    0.10994E+01    0.14614E+00   -0.13758E+01    0.26023E+02    0.29993E+02
     endif
     if ( [fittype] .eq. 'data') then
	vec/read bgw,bgc,bgerr [field]_bg_43396.dat
	vec/write  bgw,bgc,bgerr ' ' 3(f15.5,1x)
     endif
   endif
endif
* calc dilution factor ( no theta cut)
*
zone 1 2
nbin=40
wlo=.7
whi=1.100
step=([whi]-[wlo])/[nbin] 
cdir //pawc
hist/del 0
cdir //all
ttit= W  all runs [field] [name]
title_gl [ttit] 
1dhisto 1011 [ttit] [nbin]  [wlo] [whi] 
ntu/proj 1011 //all/9010.w hsp>.01.and.0.8<hsshtrk/hsp<1.4.and.abs(hsdelta)<8.
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.
hist/fit 1011([fwlo]:[fwhi]) fit_back_data.f SQ 1  fcon
endif
if ([fittype] .eq. 'yield') then
vec/create fcon(1) r 1.
sigma tot=tot*[cf]
sigma etot=etot*[cf]
vec/fit vw tot etot  fit_back_data.f ! 1  fcon
endif
if ([fittype] .eq. 'g+p2') then
vec/create fpar(6) r
vec/copy fit2 fpar
vec/create fcon(1) r 1000.
hist/fit 1011([fwlo]:[fwhi]) fit_backg.f SQ 1  fcon
endif
if ([fittype] .eq. 'p4') then
vec/create fpar(4) r
vec/copy fit2 fpar
vec/create fcon(1) r 1000.
hist/fit 1011([fwlo]:[fwhi]) fit_back.f SQ 1  fcon
endif

vec/write fcon
con=fcon(1) 
*
if ([fittype] .eq. 'yield') then
con=fcon(1)
sigma yy=bgc*[con]
endif
if ([fittype] .eq. 'data') then
con=fcon(1)
sigma yy=bgc*[con]
endif
if ([fittype] .eq. 'g+p2') then
do i=1,6
p[i]=fpar([i])
enddo
con=fcon(1)
message con = [con]
sigma yy=[con]*([p1]*exp(-0.5*(vw-[p2])*(vw-[p2])/[p3]/[p3])+[p4]+[p5]*vw+[p6]*vw*vw)
endif
if ([fittype] .eq. 'p4') then
do i=1,4
p[i]=fpar([i])
enddo
con=fcon(1)
message con = [con]
sigma yy=[con]*([p1]+[p2]*x+[p3]*x*x+[p4]*x*x*x)
endif
SET PLCI 4 
pline 40 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+bgerr*bgerr*[con]*[con]/tot/tot)
endif
if ( [fittype] .eq. 'g+p2') then
sigma edf=yy/tot/tot*etot
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] 
*
chain -all
enddo