macro angle_cdc
* 
* Plot approximation to plot momentum resolution by scaling the
*      results of momres to other scattering angles
* 02/02/07  ES
*
* set options
*
hi/create/title_global 'Scaling results of MOMRES to other angles'
option ndate
option nbox
*set stat 1111111
set stat 111
option nstat
set fit 111
option fit
option grid
zone
*
* plotting options
*
set *
*set xmgl 4.
*set ymgl 4.
*set asiz 0.4
*set xlab 2.
*set ylab 1.
*set xsiz 20.
*set xmgl 3.
*set ymgl 3.
*set ysiz 20.
*set gsiz 0.4
*
* set font definitions to bold roman 
*
set CFON -21
set GFON -21
set LFON -21
set TFON -21
set VFON -21
set txfp -21
set SMGU 0.02
set SMGR 0.02
set CSIZ 0.33
set VSIZ 0.25
set TSIZ 0.35
set YHTI 0.9
set HWID 3.0
set BWID 3.0
*
* open metafile
*
for/file 66 angle_cdc.ps
meta 66 -111
*for/file 66 angle_cdc.eps
*meta 66 -113
*
zone 2 2 
*
npts1 = 101
npts = [npts1] - 1
*
*
*
vector/create p([npts1]) R [npts1]*0.
vector/create beta([npts1]) R [npts1]*0.
vector/create delta([npts1]) R [npts1]*0.
vector/create dth([npts1]) R [npts1]*0.
vector/create dx([npts1]) R [npts1]*0.
*
sigma j=array([npts1],0#[npts1]-1)
*
pmin = 0.1
pmax = 2
mass = 0.140
alpha = 60
cosa = $sigma(cos((90-[alpha])*3.14159/180.))
mass2 = [mass]*[mass]
message 'pmin' [pmin] 'pmax' [pmax] 'mass2' [mass2]
message 'cosa' [cosa]
wait
*
dels1=0.0086
dels2=0.0111
dths1=0.00169
dths2=0.00234
dxs1=294.4
dxs2=300.3
*
sigma p = [pmin] + j*([pmax]-[pmin])/[npts]
sigma beta = p/(sqrt(p**2+[mass2]))
sigma delta = sqrt(([dels1]*p)**2 + ([dels2]/beta)**2)
sigma delta = delta/[cosa]
sigma dth = sqrt([dths1]**2 + ([dths2]/(p*beta))**2)
sigma dx = sqrt ([dxs1]**2 + ([dxs2]/(p*beta))**2)
*vector/print p
*vector/print beta
*vector/print delta
*vector/print dth
*vector/print dx
wait
*
csize = 0.05
igset chhe 0.25
*itx 0.05 0.09 'E?b! = 5.754 GeV'
*itx 0.05 0.07 'W = 2.4 GeV'
*itx 0.05 0.05 'Q^2! = 2.5 GeV^2'
* 
xmin = 0
xmax = [pmax]
ymin = 0
ymax = 0.04
set plci 2
set pmci 2
hplot/null [xmin] [xmax] [ymin] [ymax]
hplot/symbol p delta [npts1] 20 [csize]
hplot/atitle 'p (GeV/c)' '[D]p/p'
exe window#push
itx 0.5 0.9 'Theta=' [alpha]
exe window#pop
*
xmin = 0
xmax = [pmax]
ymin = 0
ymax = 0.01
hplot/null [xmin] [xmax] [ymin] [ymax]
hplot/symbol p dth [npts1] 20 [csize]
hplot/atitle 'p (GeV/c)' '[DQ] (radians)'
exe window#push
*itx 0.5 0.9 'Theta' [alpha]
itx 0.5 0.9 'Theta=90'
exe window#pop
*
set plci 2
set pmci 2
xmin = 0
xmax = [pmax]
ymin = 0
ymax = 1000
hplot/null [xmin] [xmax] [ymin] [ymax]
hplot/symbol p dx [npts1] 20 [csize]
hplot/atitle 'p (GeV/c)' '[D]x ([m]m)'
exe window#push
*itx 0.5 0.9 'Theta' [alpha]
itx 0.5 0.9 'Theta=90'
exe window#pop
wait
*
* plot 
*
vector/create theta([npts1]) R [npts1]*0.
vector/create costh([npts1]) R [npts1]*0.
vector/create delta([npts1]) R [npts1]*0.
vector/create dth([npts1]) R [npts1]*0.
vector/create dx([npts1]) R [npts1]*0.
*
* resolution at fixed momentum vs angle
*
pmom = 1.
dels1=0.0075
dels2=0.0121
dths1=0.00171
dths2=0.00263
dxs1=294.7
dxs2=349.4
*
norm = 90.
mass = 0.140
beta = $sigma([pmom]/(sqrt([pmom]**2+[mass2])))
mess ' mass' [mass] ' pmom' [pmom] ' beta' [beta]
*
zone 2 2
thmax = 150
thmin = 40
sigma theta = [thmin] + j*([thmax]-[thmin])/[npts]
sigma costh = cos((theta-[norm])*3.14159/180.)
sigma delta = sqrt(([dels1]*[pmom])**2 + ([dels2]/[beta])**2)/costh
*sigma dth = sqrt([dths1]**2 + ([dths2]/([pmom]*beta))**2)
*sigma dx = sqrt ([dxs1]**2 + ([dxs2]/([pmom]*beta))**2)
*
vector/print theta
vector/print costh
*
set plci 2
set pmci 2
xmin = [thmin]
xmax = [thmax]
ymin = 0
ymax = 0.04
*
hplot/null [xmin] [xmax] [ymin] [ymax]
hplot/symbol theta delta [npts1] 20 [csize]
hplot/atitle 'Theta (degrees)' '[D]p/p'
exe window#push
itx 0.5 0.9 'p=' [pmom]
exe window#pop
*
close 66
*
exitm 
return