macro sig_cdc
* 
* Plot approximation to plot momentum resolution 
*      as a function of position resolution in chambers
* 02/07/07  ES
*
* set options
* 
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 sig_cdc.ps
meta 66 -111
*for/file 66 sig_cdc.eps
*meta 66 -113
*
zone 2 2 
*
npts1 = 301
npts = [npts1] - 1
*
*
*
vector/create sig([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)
*
sigmin = 0
sigmax = 300
mass = 0.140
mass2 = [mass]*[mass]
message 'pmin' [pmin] 'pmax' [pmax] 'mass2' [mass2]
wait
*
dels1=0.0086
dels2=0.0111
dths1=0.00169
dths2=0.00234
dxs1=294.4
dxs2=300.3
*
sig0 = 150
*
* setup plots
*
set gsiz 0.35
set ygti 1.0
title 'MOMRES CDC 90 deg resolution vs measurement [s]'
delp = 0.4
xmin = 0
xmax = [sigmax]
ymin = 0
ymax = 0.04
hplot/null [xmin] [xmax] [ymin] [ymax]
*
* loop over p
*
do j=1,3
p = [j]*[delp]
*
beta = $sigma([p]/(sqrt([p]**2+[mass2])))
sigma sig = [sigmin] + j*([sigmax]-[sigmin])/[npts]
sigma delta = sqrt(([dels1]*[p]*sig/[sig0])**2 + ([dels2]/[beta])**2)
sigma dth = sqrt(([dths1]*sig/[sig0])**2 + ([dths2]/([p]*[beta]))**2)
sigma dx = sqrt (([dxs1]*sig/[sig0])**2 + ([dxs2]/([p]*[beta]))**2)
vector/print sig
*vector/print delta
*vector/print dth
*vector/print dx
wait
*
csize = 0.05
igset chhe 0.25
* 
set plci [j]+1
set pmci [j]+1
hplot/symbol sig delta [npts1] 20 [csize]
hplot/atitle 'Measurement [s] ([m]m)' '[D]p/p'
exe window#push
y = 1 - [j]*0.1
hplot/key 0.2 [y] 20 'p='//[p]
exe window#pop
enddo
wait
*
* setup plots
*
xmin = 0
xmax = [sigmax]
ymin = 0
ymax = 0.015
hplot/null [xmin] [xmax] [ymin] [ymax]
*
* loop over p
*
do j=1,3
p = [j]*[delp]
*
beta = $sigma([p]/(sqrt([p]**2+[mass2])))
sigma sig = [sigmin] + j*([sigmax]-[sigmin])/[npts]
sigma delta = sqrt(([dels1]*[p]*sig/[sig0])**2 + ([dels2]/[beta])**2)
sigma dth = sqrt(([dths1]*sig/[sig0])**2 + ([dths2]/([p]*[beta]))**2)
sigma dx = sqrt (([dxs1]*sig/[sig0])**2 + ([dxs2]/([p]*[beta]))**2)
*
set plci [j]+1
set pmci [j]+1
hplot/symbol sig dth [npts1] 20 [csize]
hplot/atitle 'Measurement [s] ([m]m)' '[DQ] (radians)'
exe window#push
y = 1 - [j]*0.1
hplot/key 0.2 [y] 20 'p='//[p]
exe window#pop
enddo
wait
*
* setup plots
*
xmin = 0
xmax = [sigmax]
ymin = 0
ymax = 2000
hplot/null [xmin] [xmax] [ymin] [ymax]
*
* loop over p
*
do j=1,3
p = [j]*[delp]
beta = $sigma([p]/(sqrt([p]**2+[mass2])))
sigma sig = [sigmin] + j*([sigmax]-[sigmin])/[npts]
sigma delta = sqrt(([dels1]*[p]*sig/[sig0])**2 + ([dels2]/[beta])**2)
sigma dth = sqrt(([dths1]*sig/[sig0])**2 + ([dths2]/([p]*[beta]))**2)
sigma dx = sqrt (([dxs1]*sig/[sig0])**2 + ([dxs2]/([p]*[beta]))**2)
*
set plci [j]+1
set pmci [j]+1
hplot/symbol sig dx [npts1] 20 [csize]
hplot/atitle 'Measurement [s] ([m]m)' '[D]x ([m]m)'
exe window#push
y = 1 - [j]*0.1
hplot/key 0.2 [y] 20 'p='//[p]
exe window#pop
enddo
*
* compute at mid range
*
jj = 2
p = [jj]*[delp]
beta = $sigma([p]/(sqrt([p]**2+[mass2])))
sigma sig = [sigmin] + j*([sigmax]-[sigmin])/[npts]
sigma delta = sqrt(([dels1]*[p]*sig/[sig0])**2 + ([dels2]/[beta])**2)
sigma dth = sqrt(([dths1]*sig/[sig0])**2 + ([dths2]/([p]*[beta]))**2)
sigma dx = sqrt (([dxs1]*sig/[sig0])**2 + ([dxs2]/([p]*[beta]))**2)
*
ndx1 = [npts]*100/([sigmax]-[sigmin])
ndx2 = [npts]*150/([sigmax]-[sigmin])
message 'ndx1' [ndx1] 'ndx2' [ndx2]
rdel = delta([ndx2])/delta([ndx1])
rth = dth([ndx2])/dth([ndx1])
rdx = dx([ndx2])/dx([ndx1])
message 'rdel' [rdel] 'rth' [rth] 'rdx' [rdx]
igset chhe 0.35
exe window#push
itx 1.2 0.9 'Comparisons for p =' $eval([jj]*[delp])
itx 1.2 0.8 'Ratio [s]=150[m]m/[s]=100[m]m'
itx 1.2 0.6 '[D]p/p ratio='//[rdel]
itx 1.2 0.5 '[DQ]   ratio='//[rth]
itx 1.2 0.4 '[D]x   ratio='//[rdx]
exe window#pop
*
close 66
*
exitm 
return