macro confronti val1='@%^' val2='@%^' SAMESET=0
*
* val1 = filename 1
* val2 = filename 2
*
* SAMESET >0 -> i set sono gli stessi dati/mc -> 100% correlati (es. ANI/ANII)
* SAMESET =0 -> i set sono indipendenti (es. S456/S123) 

if [val1] = '@%^' .or. [val2] = '@%^' then
 message 
 message USAGE: 'exec confronti <VAL1> <VAL2> [SAMESET]'
 message
 exitm
endif

exec leggi [val2]
* rename adds a 'r' in front of the val2 array names
exec rename
exec leggi [val1]
*
* compare takes the ratios rARRAY/ARRAY etc.
* open symbols rARRAY
* full symbols ARRAY
tit '(Full=S23GAF1-S2-7.56, Open=S456)'
exec compare [SAMESET]
*
return

macro leggi val

* normal,slope,err_normal,err_slope,chi2,n_punti_in fit expon. globale
v/cre agl(1)
v/cre bgl(1)
v/cre eagl(1)
v/cre ebgl(1)
v/cre echi2(1)
v/cre enpoin(1)

* parametri a,b,c,errori a,b,c, ch2i,#punti in fit a*exp(b*pt2+c*pt2)
v/cre a2gl(1)
v/cre b2gl(1)
v/cre c2gl(1)
v/cre ea2gl(1)
v/cre eb2gl(1)
v/cre ec2gl(1)
v/cre e2chi2(1)
v/cre e2npoin(1)

v/cre nn(2)
v/cre xl(2)
v/cre yl(2)
* dimensioni in xL,pt2; limiti di integrazione in xL e pt2
v/read NN(1),NN(2),XL(1),XL(2),YL(1),YL(2)   [val] !  o

nx=nn(1)
ny=nn(2)
v/cre xlc([nx])
v/cre dxlc([nx])
v/cre datax([nx])
v/cre sigmax([nx])
v/cre errsix([nx])
v/cre accx([nx])
v/cre eracx([nx])
v/cre recwx([nx])
v/cre gawx([nx])
v/cre fpix([nx])
v/cre erpix([nx])
v/cre ftltx([nx])
v/cre ertlx([nx])
v/cre ovlxl([nx])
v/cre ovlxle([nx])

v/cre ylc([ny])
v/cre dylc([ny])
v/cre datay([ny])
v/cre sigmay([ny])
v/cre errsiy([ny])
v/cre accy([ny])
v/cre eracy([ny])
v/cre recwy([ny])
v/cre gawy([ny])
v/cre fpiy([ny])
v/cre erpiy([ny])
v/cre ftlty([ny])
v/cre ertly([ny])
v/cre ovpt([ny])
v/cre ovpte([ny])

* dati prima di correzioni
v/cre data([nx],[ny])
v/cre acc([nx],[ny])
v/cre errac([nx],[ny])

* dati corretti
v/cre dat([nx],[ny])
v/cre erdat([nx],[ny])
v/cre sigmaxy([nx],[ny])
v/cre errxy([nx],[ny])
v/cre use2([nx],[ny])
v/cre xfit([nx])
v/cre a([nx])
v/cre aer([nx])
v/cre b([nx])
v/cre ber([nx])
v/cre chi2([nx])
v/cre npoin([nx])

* centri dei bins xL, semilarghezza, dsigma/dxL, errore, accettanza xL, errore relativo
v/read xlc,dxlc,datax,sigmax,errsix,accx,eracx [val] ! ' '

* mc rec, mc gen in xL
v/read recwx,gawx         [val] ! ' '

* 1-fraz(pi,K), errore, correz. tlt, errore, overlay XL, errore 
v/read fpix,erpix,ftltx,ertlx,ovlxl,ovlxle     [val] ! ' '

* centri dei bins pt2, semilarghezza, dsigma/dpt2, errore, accettanza pt2, errore relativo
v/read ylc,dylc,datay,sigmay,errsiy,accy,eracy [val] ! ' '

* mc rec, mc gen in pt2
v/read recwy,gawy         [val] ! ' '

* 1-fraz(pi,K), errore, correz. tlt, errore, overlay pt2, errore
v/read fpiy,erpiy,ftlty,ertly,ovpt,ovpte       [val] ! ' '


if [old] = 0 then
* parametri, errori, #punti fit singolo exp
 v/read agl,eagl,bgl,ebgl,echi2,enpoin          [val] ! ' '

* parametri, errori, #punti fit exp quadratico
 v/read a2gl,ea2gl,b2gl,eb2gl,c2gl,ec2gl,e2chi2,e2npoin [val] ! ' '
endif

* dati prima di correzioni e sottrazioni di bins
v/read data     [val] ! ' '

* accettanza (prima di sottrazione)
v/read acc      [val] ! ' '

* errore accettanza
v/read errac    [val] ! ' '

* dati dopo correzioni e sottrazione di bins
v/read dat      [val] ! ' '
* errore dati
v/read erdat    [val] ! ' '

* sezione 2dif, dopo sottrzione dei bins
v/read sigmaxy  [val] ! ' '
* errore sigmaxy
v/read errxy    [val] ! ' '

* matrice sottrazione bins nei bins dell'analisi
v/read use2     [val] ! ' '

do i=1,[nx]
 i1=[i]+1
*         xL      a      da       b      db  --> fit a*exp(b*pt2)
 v/read xfit([i]),a([i]),aer([i]),b([i]),ber([i]),chi2([i]),npoin([i]) [val] ! ' '
enddo

v/cre eof(1)
v/read eof [val] ! 'c'
end_of_file=eof(1)
if [end_of_file] .ne. 99999 then
 mess *** ERROR, END_OF_FILE EXPECTED!!!!
endif
v/del eof

return

macro butta_fuori ofile
fmt=(1x,i5,4F10.3,1x,2i5)
nx=$VDIM(xlc,1)
ny=$VDIM(ylc,1)
nxy=[nx]*[ny]
v/cre out(4) r
v/cre use(1) i
v/cre ind(2) i
do i=1,[nx]
 v/inp ind(1) [i]
 do j=1,[ny]
   v/inp ind(2) [j]
   ij=[i]*[j]
   sxy=sigmaxy([i],[j])
   exy=errxy([i],[j])
   if [sxy] > 0 then
    v/inp use(1) 1
   else
    v/inp use(1) -1
   endif
   v/inp out(1) xlc([i])
   v/inp out(2) ylc([j])
   v/inp out(3) [sxy]
   v/inp out(4) [exy]
   if [ij] = 1 then
    v/write use,out(1),out(2),out(3),out(4),ind(1),ind(2) [ofile] [fmt] O
   elseif [ij] = [nxy] then
    v/write use,out(1),out(2),out(3),out(4),ind(1),ind(2) [ofile] [fmt] C
   else
    v/write use,out(1),out(2),out(3),out(4),ind(1),ind(2) [ofile] [fmt] ' '
   endif
  enddo
enddo
return



macro rename

* normal,slope,err_normal,err_slope,chi2,n_punti_in fit expon. globale
v/copy agl ragl
v/copy bgl rbgl
v/copy eagl reagl
v/copy ebgl rebgl
v/copy echi2 rechi2
v/copy enpoin renpoin

v/copy a2gl ra2gl
v/copy b2gl rb2gl
v/copy c2gl rc2gl
v/copy ea2gl rea2gl
v/copy eb2gl reb2hl
v/copy ec2gl rec2gl
v/copy e2chi2 re2chi2
v/copy e2npoin re2npoin

v/copy xlc rxlc
v/copy dxlc rdxlc
v/copy datax rdatax
v/copy sigmax rsigmax
v/copy errsix rerrsix
v/copy accx   raccx
v/copy eracx  reracx
v/copy recwx  rrecwx
v/copy gawx   rgawx
v/copy fpix   rfpix
v/copy erpix  rerpix
v/copy ftltx  rftltx
v/copy ertlx  rertlx
v/copy ovlxl  rovlxl
v/copy ovlxle rovlxle

v/copy ylc    rylc
v/copy dylc   rdylc
v/copy datay  rdatay
v/copy sigmay rsigmay
v/copy errsiy rerrsiy
v/copy accy   raccy
v/copy eracy  reracy
v/copy recwy  rrecwy
v/copy gawy   rgawy
v/copy fpiy   rfpiy
v/copy erpiy  rerpiy
v/copy ftlty  rftlty
v/copy ertly  rertly
v/copy ovpt   rovpt
v/copy ovpte  rovpte

* dati prima di correzioni
v/copy data  rdata
v/copy acc   racc
v/copy errac rerrac

* dati corretti
v/copy dat   rdat
v/copy erdat rerdat
v/copy sigmaxy rsigmaxy
v/copy errxy   rerrxy
v/copy use2    ruse2
v/copy xfit    rxfit
v/copy a ra
v/copy aer raer
v/copy b rb
v/copy ber rber
v/copy chi2 rchi2
v/copy npoin rnpoin

return

macro compare sameset=0 
n=$VDIM(xlc,1)
m=$VDIM(ylc,1)
opt pto
zone 1 2
null 0.18 1.1 0. 500.
hplot/err xlc sigmax dxlc errsix [n] 20 0.2
hplot/err xlc rsigmax dxlc rerrsix [n] 24 0.3
*exec setkey 0.1 0.90  24 ! LR
*exec setkey 0.1 0.80 20 ! GB
atit 'x?L!' 'd[s]/dx?L! (nb)'
sigma r=rsigmax/sigmax
if [sameset] > 0 then
 sigma dr=rerrsix*r*sqrt(1/sigmax**2+1/rsigmax**2-2/(sigmax*rsigmax))
else
 sigma dr=r*sqrt((rerrsix/rsigmax)**2+(errsix/sigmax)**2)
endif
* null 0.18 1.1 0.8 1.2
null 0.18 1.1 0.1 2.
hplot/err xlc r dxlc dr [n] 20 0.15
dline 0.18 1.1 1. 1.
atit 'x?L'
* exec setcom 0.2 0.8 'RATIO LR/GB' 0.4
null 0 0.5 0. 500.
hplot/err ylc sigmay  dylc  errsiy [m] 20 0.2
hplot/err ylc rsigmay dylc rerrsiy [m] 24 0.3
*exec setkey 0.9 0.90  24 ! LR
*exec setkey 0.9 0.80 20 ! GB
atit 'P?T!!^2!' 'd[s]/dP?T!!^2! (nb/GeV^2!)'
sigma r=rsigmay/sigmay
if [sameset] > 0 then
 sigma dr=rerrsiy*r*sqrt(1/sigmay**2+1/rsigmay**2-2/(sigmay*rsigmay))
else
 sigma dr=r*sqrt((rerrsiy/rsigmay)**2+(errsiy/sigmay)**2)
endif
null 0. 0.5 0.1 2.0
hplot/err ylc r dylc dr 15 20 0.15
dline 0. 0.5  1. 1.
atit 'P?T!!^2!' 
*
* exec setcom 0.2 0.8 'RATIO LR/GB' 0.4
*
zone 2 2
v/cop RSIGMAXY DUM
sigma RATIO=DUM/SIGMAXY

if [SAMESET] > 0 then
 sigma ERR=RATIO*RERRXY*SQRT(1/SIGMAXY**2+1/DUM**2-2/(SIGMAXY*DUM))
else
 sigma ERR=RATIO*SQRT((ERRXY/SIGMAXY)**2+(RERRXY/DUM)**2)
endif

v/cre R([n])
v/cre DR([n])
v/cop R A
v/cop R DA
v/cop R B
v/cop R DB
sigma XA=VMIN(XLC-DXLC)*0.9
sigma XB=VMAX(XLC+DXLC)*1.1
x1=XA(1)
x2=XB(1)
y1=0.
y2=800.
rmin=0.
rmax=2.
do j=1,[m]
 do i=1,[n]
   v/inp R([i])   RATIO([i],[j])
   v/inp DR([i])  ERR([i],[j])
   v/inp A([i])   RSIGMAXY([i],[j])
   v/inp B([i])   SIGMAXY([i],[j])
   v/inp DA([i])   ERRXY([i],[j])
   v/inp DB([i])   RERRXY([i],[j])
  enddo
  null [x1] [x2] [y1] [y2]
  hplot/err XLC A DXLC DA [n] 20 0.1
  hplot/err XLC B DXLC DB [n] 24 0.2  
  atit 'x?L!'
  null [x1] [x2] [rmin] [rmax]
  hplot/err XLC R DXLC DR [n] 20 0.2
*  dline [x1] [x2] 1. 1.
  dline [x1] [x2] 1.3 1.3
  dline [x1] [x2] 0.7 0.7
  atit 'x?L!'
enddo

return