MACRO PAWTEST2
*
*   macro derived from an initial macro by Rob Veenhof
*
histogram/file 21 pawdemo.hbook
histogram/delete 0
Set *
Option *
Size 20 24
Set GSIZ 1.0
Set YGTI 0.6
set *fon -40
Title 'PAW test number 2'
Option FIT
*Set BWID 8
*Set PWID 8
Zone 1 2
axes='log'
lab=1
colour=1
*-----------------------------------------------------------------------
* Obtain the counts of the various resonances.
*-----------------------------------------------------------------------
exec rescnt
*-----------------------------------------------------------------------
* Create an appropriate directory structure.
*-----------------------------------------------------------------------
cd //PAWC
mdir BACK
cd //PAWC/BACK
1d 10001 'mix'               150 0 1.5
1d 20001 'mix err2'          150 0 1.5
1d 10002 'eta Dalitz'        150 0 1.5
1d 20002 'eta Dalitz err2'   150 0 1.5
1d 10003 'etap Dalitz'       150 0 1.5
1d 20003 'etap Dalitz err2'  150 0 1.5
1d 10004 'r/w'               150 0 1.5
1d 20004 'r/w err2'          150 0 1.5
1d 10005 'phi'               150 0 1.5
1d 20005 'phi err2'          150 0 1.5
1d 10006 'omega pi0'         150 0 1.5
1d 20006 'omega pi0 err2'    150 0 1.5
1d 10007 'eta mumu'          150 0 1.5
1d 20007 'eta mumu err2'     150 0 1.5
1d 10008 'rho pipi'          150 0 1.5
1d 20008 'rho pipi err2'     150 0 1.5
1d 10009 'phi K K'           150 0 1.5
1d 20009 'phi K K  err2'     150 0 1.5
1d 10010 'charm'             150 0 1.5
1d 20010 'charm err2'        150 0 1.5
1d 10011 'BS'                100 0 2.0
1d 20011 'BS err2'           100 0 2.0
*-----------------------------------------------------------------------
* Normalise the backgrounds.
*-----------------------------------------------------------------------
exec mixcal2
exec etacal
exec etapcal
exec veccal
exec hadcal
exec charm
exec brems
*-----------------------------------------------------------------------
* Plot the results.
*-----------------------------------------------------------------------
cd //pawc/back

vec/cr a(150)
add 20001 20002 20000
add 20000 20003 20000
add 20000 20004 20000
add 20000 20005 20000
add 20000 20006 20000
add 20000 20007 20000
add 20000 20008 20000
add 20000 20009 20000
add 20000 20010 20000
h/get/cont 20000 a
sigma a=sqrt(a)
h/put/cont 20000 a
vec/del a

add 10001 10002 10000
add 10000 10003 10000
add 10000 10004 10000
add 10000 10005 10000
add 10000 10006 10000
add 10000 10007 10000
add 10000 10008 10000
add 10000 10009 10000
add 10000 10010 10000

set xlab 1.3
if [colour]=1 then
     set fwid 3
     set hwid 3
     set bwid 3
endif
sub 10000 20000 30001
add 10000 20000 30002
add 20000 20000 30000
v/cr temp(150)
get/cont 30001 temp ; v/inp temp(1:20) 20*0
put/cont 30001 temp
max 30001 700
set hwid 4
if [axes]='log' then
     opt logy
     set xval 0.4
     min 30001 1
else
     opt liny
     set xval 0.1
endif
h/pl 30001 K
atitle 'm?[mm]! "M#GeV"N#' 'Events / 20 MeV'
set dmod 1
set hcol 1005
get/cont 30000 temp ; v/inp temp(1:21) 21*0
put/cont 30000 temp
set dmod 1
h/pl 30000         +
set hcol 0
set dmod 1
h/pl 30001      s
set dmod 1
h/pl 30002(20:) s
goto skip2
skip:
set dmod 1
h/pl 10000(20:) s
skip2:

If [lab]=1 Then
     vec/cr par(10)
     if [colour]=1 then
          igset txci 7; set fcol 7
     endif
     set dmod 1
     set fwid 8
     h/fit 10001(25:150) p7   '0q'        ; h/pl 10001(23:150) 'func,s'
     text 0.55 130 'hadrons'        0.3 -22 l
*
     if [colour]=1 then
          igset txci 2; set fcol 2
     endif
     set dmod 1
     h/fit 10002(23:55)  p7   '0q'        ; h/pl 10002(23:150) 'func,s'
     text 0.31  145 '[c]"5#[mmg]'    0.3   -46 c
*
     if [colour]=1 then
          igset txci 2; set fcol 2
     endif
     set dmod 1
     h/fit 10003(23:87)  p7   '0q'        ; h/pl 10003(23:150) 'func,s'
     text 0.27  35 '[c]^,!"5#[mmg]' 0.3 -27 l
*
     if [colour]=1 then
          igset txci 3; set fcol 3
     endif
     set dmod 1
     smooth 10004 q1n;             smooth 10004  qn
      h/pl 10004(23:150) 'func,s'
     text 0.81  90 '[r/w]"5#[mm]'   0.3  90 c
*
     if [colour]=1 then
          igset txci 3; set fcol 3
     endif
     set dmod 1
     h/fit 10005(90:120) g    '0qw'       ; h/pl 10005(23:150) 'func,s'
     text 1.06   6 '[f]"5#[mm]'     0.3  90 c
*
     if [colour]=1 then
          igset txci 6; set fcol 6
     endif
     set dmod 1
     h/fit 10006(23:70)  p7   '0ql'       ; h/pl 10006(23:150) 'func,s'
     text 0.53  77 '[w]"5#[mmp]^0!' 0.3   0 c
*
     if [colour]=1 then
          igset txci 3; set fcol 3
     endif
     set dmod 1
     h/fit 10007(40:70)  g    '0qw'       ; h/pl 10007(23:150) 'func,s'
     text 0.58 2.5 '[c]"5#[mm]'     0.3  90 c
*
     if [colour]=1 then
          igset txci 4; set fcol 4
     endif
     vec/inp par(1:7) 1 0.7 0.1 0 0 0 0
     set dmod 1
     h/fit 10008(23:95)  g+p3 '0qw' 7 par ; h/pl 10008(23:150) 'func,s'
     text 0.73 6.7   '[r]"5#[pp]'     0.3 -55 c
*
     vec/inp par(1:7) 1 0.7 0.1 0 0 0 0
     if [colour]=1 then
          igset txci 4; set fcol 4
     endif
     set dmod 1
     h/fit 10009(23:60)  g+p3 '0qw' 7 par ; h/pl 10009(23:150) 'func,s'
     text 0.41 1.8 '[f]"5#[KK]'     0.3  -50 c
*
     if [colour]=1 then
          igset txci 1; set fcol 1
     endif
     set dmod 1
     smooth 10010 qn;          h/pl 10010(23:150) 'func,s'
     text 1.24 1.8 'charm'          0.3  -19 c
*
     if [colour]=1 then
          igset txci 1; set fcol 1
     endif
     set dmod 1
     h/fit 10011(14:60) e     '0qw'       ; h/pl 10011(14:60) 'func,s'
     text 0.30 10   'Bremsstrahlung' 0.3 -55 c
*
     vec/del par
     if [colour]=1 then
          igset txci 1; set fcol 0; set fwid 1; igset pass 1
     endif
Else
     if [colour]=1 then
          igset txci 7; set hcol 7
          set hwid 3
     endif
     set dmod 1
     h/pl 10001(23:150) s
     if [colour]=1 then
          igset txci 2; set hcol 2
     endif
     set dmod 1
     h/pl 10002(23:150) s
     if [colour]=1 then
          igset txci 2; set hcol 2
     endif
     set dmod 1
     h/pl 10003(23:150) s
     if [colour]=1 then
          igset txci 3; set hcol 3
     endif
     set dmod 1
     h/pl 10004(23:150) s
     if [colour]=1 then
          igset txci 3; set hcol 3
     endif
     set dmod 1
     h/pl 10005(23:150) s
     if [colour]=1 then
          igset txci 6; set hcol 6
     endif
     set dmod 1
     h/pl 10006(23:150) s
     if [colour]=1 then
          igset txci 3; set hcol 3
     endif
     set dmod 1
     h/pl 10007(23:150) s
     if [colour]=1 then
          igset txci 4; set hcol 4
     endif
     set dmod 1
     h/pl 10008(23:150) s
     if [colour]=1 then
          igset txci 4; set hcol 4
     endif
     set dmod 1
     h/pl 10009(23:150) s
     if [colour]=1 then
          igset txci 1; set hcol 1
     endif
     set dmod 1
     h/pl 10010(23:150) s
     if [colour]=1 then
          igset txci 1; set hcol 1
     endif
     set dmod 1
     h/pl 10011         s
     if [colour]=1 then
          igset txci 1; set hcol 0
          set hwid 1
     endif
Endif

cd //lun21/mm_mmg
vec/cr xma(75); vec/cr xme(75); vec/cr yma(75); vec/cr yme(75)
rebin 2 xma(1:75) yma(1:75) xme(1:75) yme(1:75) 75 1 150
sigma yma=2*yma
sigma yme=2*yme
set lwid 2
hpl/errors xma(11:75) yma(11:75) xme(11:75) yme(11:75)  65 20  0.25

close 21
v/del *
set *
opt *
zon 1
title
cd //pawc/back; h/del 0; cd //pawc; h/del 0; ddir BACK
leave:
return

MACRO RESCNT
cd //lun21/mm_mmg
exec vector
exec pseudo
cd //pawc
neta=neta
nrw=nrw
nphi=nphi
message 'Resonance count: eta=' [neta] 'r/w=' [nrw] 'phi=' [nphi]
h/del 0
return

MACRO VECTOR
vec/cr a(150)
vec/cr f(150)
h/get/cont 2 a
sigma b=a
vec/inp a(50:60) 11*0
vec/inp a(70:86) 17*0
vec/inp a(92:112) 21*0
h/put/cont 2 a
h/fit 2(35:140) p4 'w0q'
h/get/func 2 f
sigma a=b-f
h/put/cont 2 a
vec/cr par(3)
vec/cr dum(3)
vec/cr epar(3)
h/fit 2(55:67) g '0q' 0 par dum dum epar
message 'Measured eta -> mm mass:'; sigma print par(2)
sigma neta=par(1)*par(3)*sqrt(2*pi)/0.010
message 'Events in eta -> mu mu: '; sigma print neta
h/fit 2(70:93) g '0q' 0 par dum dum epar
message 'Measured r/w -> mm mass:'; sigma print par(2)
sigma nrw=par(1)*par(3)*sqrt(2*pi)/0.010
h/fit 2(100:112) g '0q' 0 par dum dum epar
message 'Measured phi -> mm mass:'; sigma print par(2)
sigma nphi=par(1)*par(3)*sqrt(2*pi)/0.010
sigma nrw=nrw*1.2
sigma nphi=nphi*1.2
vec/del a; vec/del f; vec/del par; vec/del dum; vec/del epar
*set bcol 1107
set hcol 1105
set fcol 4
set fwid 8
set cshi 0.02
h/plot 2
set bcol 0
set hcol 0
set fwid 1
set fcol 1
*  First zone completed
return

MACRO PSEUDO
vec/cr a(100)
vec/cr f(100)
hist=10
h/get/cont [hist] a
sigma b=a
vec/inp a(20:40) 21*0
h/put/cont [hist] a
h/fit [hist](10:70) p4 '0q'
h/get/func [hist] f
sigma a=b-f
h/put/cont [hist] a
vec/cr par(3)
vec/cr dum(3)
vec/cr epar(3)
h/fit [hist](10:50) g '0q' 0 par dum dum epar
sigma sigeta=par(2)
sigeta=sigeta
message 'Experimental eta width:' [sigeta] 'GeV'
sigma neta=par(1)*par(3)*sqrt(2*pi)/0.020
vec/del a; vec/del f; vec/del par; vec/del dum; vec/del epar
return

MACRO ETACAL
vec/cr par(3); vec/cr epar(3); vec/cr dum(3)
vec/cr a(150); vec/cr err(150); vec/cr gen(150)

cd //lun21/etaf2
h/fit 4(10:50) g '0q' 0 par dum dum epar
sigma sigeta=par(2)
sigeta=sigeta
message 'Monte-Carlo  eta width:' [sigeta] 'GeV'
sigma netamc=par(1)*par(3)*sqrt(2*pi)/0.020
h/get/cont 103 a
h/get/cont 109 gen
sigma a=a*neta/netamc
sigma err=(a*0.25)**2
sigma geneta=vsum(gen)*neta/netamc

cd //pawc/back
h/put/cont 10002 a
h/put/cont 20002 err

cd //pawc
h/del 0
vec/del par; vec/del epar; vec/del dum; vec/del a; vec/del err; vec/del gen
return

MACRO ETAPCAL
vec/cr a(150); vec/cr err(150); vec/cr gen(150)
cd //lun21/etapfv1
h/get/cont 1109 gen
h/get/cont 1103 a
sigma a=0.145*a*geneta/vsum(gen)
sigma err=(0.55*a)**2
cd //pawc/back
h/put/cont 10003 a
h/put/cont 20003 err
vec/del a; vec/del err; vec/del gen
return

MACRO VECCAL
vec/cr par(3); vec/cr dum(3); vec/cr epar(3); vec/cr a(150); vec/cr err(150)
vec/cr gen(150)

h/get/cont //lun21/ro/103 a
sigma a=2*a*nrw/vsum(a)
sigma err=(0.3*a)**2
cd //pawc/back
h/put/cont 10004 a
h/put/cont 20004 err

h/get/cont //lun21/phi/103 a
h/get/cont //lun21/phi/101 gen
sigma genphi=vsum(gen)*nphi/vsum(a)
sigma a=a*2*nphi/vsum(a)
sigma err=(0.07*a)**2
cd //pawc/back
h/put/cont 10005 a
h/put/cont 20005 err

h/get/cont //lun21/omegaf/103 a
h/get/cont //lun21/omegaf/109 gen
sigma a=a*7.69*genphi/vsum(gen)
sigma err=(0.45*a)**2
cd //pawc/back
h/put/cont 10006 a
h/put/cont 20006 err

h/get/cont //lun21/etamumu/103 a
h/get/cont //lun21/etamumu/101 gen
sigma a=a*1.35*genphi/vsum(gen)
sigma err=(1.00*a)**2
cd //pawc/back
h/put/cont 10007 a
h/put/cont 20007 err

vec/del par; vec/del dum; vec/del epar; vec/del a; vec/del err; vec/del gen
return

MACRO HADCAL
vec/cr a(150); vec/cr err(150); vec/cr gen(150)

h/get/cont //lun21/pion/103 a
h/get/cont //lun21/pion/101 gen
sigma a=a*194000*genphi/vsum(gen)
sigma err=(0.30*a)**2
cd //pawc/back
h/put/cont 10008 a
h/put/cont 20008 err

h/get/cont //lun21/kaon/103 a
h/get/cont //lun21/kaon/101 gen
sigma a=a*1830*genphi/vsum(gen)
sigma err=(0.25*a)**2
cd //pawc/back
h/put/cont 10009 a
h/put/cont 20009 err

vec/del a; vec/del err; vec/del gen
return

MACRO CHARM
cd //lun21/charm
vec/cr par(5); vec/cr c(150); vec/cr v(150)
h/fit 7103 p4 '0q' 5 par
cd //pawc
sigma c=array(150,0#1.5)
sigma v=par(1)+par(2)*c+par(3)*c**2+par(4)*c**3+par(5)*c**4
cd //pawc/back
h/put/cont 10010 v
sigma err=(0.4*v)**2
h/put/cont 20010 err
vec/del par; vec/del c; vec/del v; vec/del err
return

MACRO BREMS
vec/cr c(100)
cd //lun21/bs
h/get/cont 200 c
sigma cc=2*c*2.3e7/521
cd //pawc/back
h/put/cont 10011 cc
sigma err=(0.5*cc)**2
h/put/cont 20011 err
vec/del c; vec/del cc; vec/del err;
return

MACRO MIXCAL
cd //lun21/mm_mmg
vec/cr a(150); vec/cr pp(150); vec/cr mm(150); vec/cr err(150)
sigma ord=array(150,0.010#1.490)
sigma corr=0.524+4.071*ord-6.158*ord**2+3.640*ord**3-0.742*ord**4
h/get/cont 5 a
h/get/cont 1 pp
h/get/cont 3 mm
cd //pawc/back
sigma a=a*4*corr*sqrt(vsum(pp)*vsum(mm))/vsum(a)
sigma err=(0.1*a)**2
h/put/cont 10001 a
h/put/cont 20001 err
cd //pawc
h/del 0
vec/del a; vec/del pp; vec/del mm; vec/del err; vec/del ord; vec/del corr
return

MACRO MIXCAL2
cd //lun21/mm_mmg
vec/cr pp(150)
vec/cr mm(150)
vec/cr pm(150)
vec/cr err(150)
h/get/cont 1 pp
h/get/cont 3 mm
message calling Application SIGMA
Application SIGMA quit
   ord=array(150,0.010#1.490)
   corr=0.524+4.071*ord-6.158*ord**2+3.640*ord**3-0.742*ord**4
   pm=sqrt(pp*mm)
   pm=pm*4*corr*sqrt(vsum(pp)*vsum(mm))/vsum(pm)
   err=(0.1*pm)**2
QUIT
1d 6 ' ' 150 0 1.5
h/put/cont 6 pm
vec/cr par(8)
vec/inp par 3.3 -5 -8 45 27 -17 0 0
h/fit 6(18:150) e*p3 '0q' 7 par
h/get/fun 6 pp
vec/inp pp(1:21) 21*0
h/put/con 10001 pp
h/put/con 20001 err
cd //pawc
h/del 0
vec/del pp; vec/del pm; vec/del mm; vec/del par; vec/del err
vec/del ord; vec/del corr
return