*************************************************************
*	kumac modificata il 9/07/2009
* e' stato cambiato il grafico ophase.eps che comunque non e' corretto
*************************************************************
macro proxie file='/home/mcontalb/CLAS/newRICH/OUTPUT/0000_proxie.rz'
 close 1

 if ($FEXIST([file]) .EQ. 1) then
   h/fil 1 [file]
   n/pri 10
 else
   message ERROR: file [file] does not exist
   stopm
 endif

return
 

* **********************
* Load detector geometry 
* 
* return the vector vddg
* return global surdet = detector surface in m2 if input is in cm2
 
 macro loaddg ddfile='0000_ddg.txt'
   
   m/gl/import surdet 0.
   
   v/del vddg
   
   IF ($FEXIST([ddfile]) .EQ. 0) THEN
     MESSAGE WARNING [ddfile] detector definition file does not exist
     MESSAGE use defaults values
     V/CREA vddg(5) R 57.82 229.54 -114.77 -30. 90.
   ELSE
     message [ddfile] detector definition file loaded
     v/read vddg [ddfile]
   ENDIF
   v/pri vddg
* area (surface)
      
   da=$SIGMA((vddg(5)-vddg(4))/180.*3.14)
   s1=$SIGMA(vddg(1)*vddg(1)*[da]/2.)
   s2=$SIGMA(vddg(2)*vddg(2)*[da]/2.)
   
   surdet=([s2]-[s1])/10000.
   
 return
 
* **********************
* Load radiator geometry 
* 
* return the vector vddg
 
 macro loardg rdfile='0000_rdg.txt'
   
   m/gl/import surrad 0.
   
   v/del vrdg
   
   IF ($FEXIST([rdfile]) .EQ. 0) THEN
     MESSAGE WARNING: [rdfile] radiator definition file does not exist
     MESSAGE use defaults values
     V/CREA vrdg(4) R 50.72 201.33 -100.67 27.1
   ELSE
     message [rdfile] radiator definition file loaded
     v/read vrdg [rdfile]
   ENDIF
   v/pri vrdg
   
* area in m2 if input is cm
   
   da=$SIGMA(2.*vrdg(4)/180.*3.14)
   s1=$SIGMA(vrdg(1)*vrdg(1)*[da]/2.)
   s2=$SIGMA(vrdg(2)*vrdg(2)*[da]/2.)
   
   surrad = ([s2] - [s1])/10000.
   
 return
   
* ********************************************************
* Reconstruct the Angle in clas12 proximity rich
* 
* .... better to call lrclas
* 
* error = 1: chromatic
* error = 2: chromatic + emission
* error = 3: chromatic + emission + localization
 
 macro rclas wf=1.5 ef=0.6 outf='' p_gev0=2.0 p_gev1=6.0 error=3 nbin=40 ddfile='0000_ddg.txt' outps='cher_vs_mom.eps' giro=1

   opt nbox
   set *fon 60
   set chhe .4
   set csiz .5
   
   set gsiz .5
   set ygti 1.0

   set asiz .45
   set vsiz .4
   set xval .2

   opt utit
   
   message '==========================================='
   message 'FREON THICKNESS = ' [wf] 'cm'
   message 'DETECTOR GLOBAL EFFICIENCY = ' [ef]
   message '==========================================='

   xmas1=0.93828
   xmas2=0.493677
   xmas3=0.13957

   h/del *
   
   if ($VEXIST(vddg) .EQ. 0) then
     exec proxie#loaddg [ddfile]
   endif
   
* histogram min, max angle

   m/gl/import ns1 0.
   m/gl/import ns2 0.
   m/gl/import ns3 0.

   m/gl/import sk1 0.
   m/gl/import sk2 0.
   m/gl/import sk3 0.

   m/gl/import spi1 0.
   m/gl/import spi2 0.
   m/gl/import spi3 0.
   
   m/gl/import mk3 0.
   m/gl/import mpi3 0.
   
   title 'Cherenkov Angles vs Momentum'
   zone 1 4
   opt nfit
   opt nsta
   set chhe .35
   set fwid 3.
   set plci 2
   set ywin 1.
   set ylab .6
   set yval .1
   p0=.53
   p1=.57
   k0=.64
   k1=.66
   pi0=.67
   pi1=.75

   par1=14 
   par2=11 
   par3=8  

   
* n/plot 10.theta0g 25.<theta0g<55.
   
   cut $1 ([p_gev0]<sqrt(px*px+py*py+pz*pz)<[p_gev1])
   
   ymin=0.5
   ymax=0.725
   yym=0.005

   message 'limits' [p_gev0] [p_gev1] [ymin] [ymax]
   
   h/cre/profile 31 'p Cherenkov vs Momentum' [nbin] [p_gev0] [p_gev1] [ymin] [ymax] S
   h/cre/profile 32 'K Cherenkov vs Momentum' [nbin] [p_gev0] [p_gev1] [ymin] [ymax] S
   h/cre/profile 33 '[p] Cherenkov vs Momentum' [nbin] [p_gev0] [p_gev1] [ymin] [ymax] S
   h/cre/profile 41 'p Cherenkov vs Momentum' [nbin] [p_gev0] [p_gev1] -[yym] [yym] S
   h/cre/profile 42 'K Cherenkov vs Momentum' [nbin] [p_gev0] [p_gev1] -[yym] [yym] S
   h/cre/profile 43 '[p] Cherenkov vs Momentum' [nbin] [p_gev0] [p_gev1] -[yym] [yym] S

   
   close 10
   fort/file 10 [outps]
   meta 10 -113
   
   null [p_gev0] [p_gev1] [ymin] [ymax]
   
* generate mask file
   
   mask/file stmask N
   
   do i=1,3
     message 'ipar' par[i] [par[i]]
     n/loop 10 ipar=[par[i]]>>stmask([i])
   enddo
   n/loop 10 prosel.f([error].,[wf],[ef])>0.&&$1>>stmask(4)
   n/loop 10 prosel.f(4.,[wf],[ef])>0.&&$1>>stmask(5)

   
* xxxx
   
   do i=1,3 
     hj=30+[i]
     hjj=50+[i]
     hjjj=40+[i]
     mj=20+2*([i]-1)
     set hcol [i]
     set hwid 5
     set pmci [i]
     set mtyp [mj]
     n/plot 10.(prosel.f([error].,[wf],[ef]))%sqrt(px*px+py*py+pz*pz) prosel.f(4.,[wf],[ef])>10&&stmask(4)&&stmask([i])  option=sprofs idh=[hj]

     v/cre v[i]([nbin]) r
     v/cre ev[i]([nbin]) r
     v/cre evv[i]([nbin]) r
     v/crea xab[i]([nbin]) r
     v/crea nv[i]([nbin]) r
     
     h/get/cont [hj] v[i]
     h/get/absc [hj] xab[i]
     h/get/error [hj] evv[i]
     
     sigma vmedang=v[i]
     n/plot 10.(prosel.f([error].,[wf],[ef]))-medang.f([i],[nbin],[p_gev0],[p_gev1])%sqrt(px*px+py*py+pz*pz) prosel.f(4.,[wf],[ef])>10&&stmask(4)&&stmask([i]) option=sprofs idh=[hjjj]
     h/get/error [hjjj] ev[i]

     n/pro [hjj] 10.sqrt(px*px+py*py+pz*pz) prosel.f(4.,[wf],[ef])>10&&stmask(4)&&stmask([i]) 
     h/get/cont [hjj] nv[i]

* - valore medio errore
     m[i]=$SIGMA(VSUM(ev[i])/NCO(ev[i]))*1000.
   enddo  

*  close 90
   
   v/crea vmom([nbin]) r   
   h/get/abscissa [hj] vmom
   
   atitle 'p (GeV/c)' '[q] (rad)'
   
   sigma nskp=2.*(v2-v1)/(ev2+ev1)
   sigma nskpi=2.*(v3-v2)/(ev3+ev2)

   message [m1] [m2] [m3]
   
*  ----------------------------------------
*  -  distribuzione errori
*  ----------------------------------------
   null 0 10 0 10 sab
   titolo=Proton 
* [s]=$UNQUOTE($FORMAT([m1],F5.2)) mr
   set pmci 1
   key 7 2. 20 [titolo] 0.5
   titolo=Kaon 
* [s]=$UNQUOTE($FORMAT([m2],F5.2)) mr
   set pmci 2
   key 7 3.5 22 [titolo] 0.5
   titolo=Pion 
* [s]=$UNQUOTE($FORMAT([m3],F5.2)) mr
   set pmci 3
   key 7 5 24 [titolo] 0.5
   
*  error 
   op=paw
   do i=1,3
     mj=20+2*([i]-1)
     set hcol [i]
     set hwid 5
     set pmci [i]
     set mtyp [mj]
     sigma evm=ev[i]*1000.
     graph [nbin] xab[i] evm [op]
     op=p
   enddo
   atitle 'p (GeV/c)' '[s]?[q]! (mrad)'
   
   
*  ----------------------------------------
*  -  Numero foto elettroni
*  ----------------------------------------
   sopt='profs'
   do i=1,3
     hj=30+[i]
     mj=20+2*([i]-1)
     set hcol [i]
     set hwid 5
     set pmci [i]
     set mtyp [mj]
* n/plot 10.(prosel.f(4.,[wf],[ef]))%sqrt(px*px+py*py+pz*pz) prosel.f(4.,[wf],[ef])>0.&&$1&&ipar=[par[i]] option=[sopt]
     n/plot 10.(prosel.f(4.,[wf],[ef]))%sqrt(px*px+py*py+pz*pz) stmask(5)&&stmask([i]) option=[sopt]
     sopt=s[sopt]
   enddo
   
   atitle 'p (GeV/c)' 'N?PE!'
   
   sopt='P'
   do i=1,3
     hj=30+[i]
     mj=20+2*([i]-1)
     set hcol [i]
     set hwid 5
     set pmci [i]
     set mtyp [mj]
* n/plot 10.(prosel.f(4.,[wf],[ef])) prosel.f(4.,[wf],[ef])>0.&&$1&&ipar=[par[i]] option=[sopt]
     n/plot 10.(prosel.f(4.,[wf],[ef])) stmask(5)&&stmask([i]) option=[sopt]
     sopt=s[sopt]
     npe[i]=$HINFO(1000000,'MEAN')
   enddo
   
   atitle 'Npe'
   
   v/del vnne
   
   v/crea vnne(3) R [npe1] [npe2] [npe3]
   
   message ' $$$ NPE $$$ ' [npe1] [npe2] [npe3]
   
   mask/close stmask
  
   close 10
   
 return
 
* ********************
* Loop on rclas
*  
* 
 
 macro lrclas file='none' nloop=1 wf=1.5 ef=0.6 outf='' p_gev0=2.0 p_gev1=6.0 error=3 nbin=40 outps='lrclas.eps' ddfile='none' rdfile='none'
   
   m/gl/create surdet 0.
   m/gl/create surrad 0.
   
   case [file] in
     (*)     message 'passo per vai'
             exec vai [file] 
     (diretto) message 'proxie diretto'
               exec proxie [file].rz
   endcase
   
   exec proxie#loaddg [ddfile]
   
   exec proxie#loardg [rdfile]

   v/crea ns1([nbin]) R
   v/crea ns2([nbin]) R
   v/crea ens1([nbin]) R
   v/crea ens2([nbin]) R
   
   msig1=0.
   msig2=0.
   msig3=0.
   
   oofile=[outps]_r.eps
   
   do k=1,[nloop]
     message '################ Lancio rclas: ' [k]
     exec proxie#rclas wf=[wf] ef=[ef] p_gev0=[p_gev0] p_gev1=[p_gev1] error=[error] nbin=[nbin] outps=[oofile] giro=[k]
     
     do i=1,3
       msig[i]=[msig[i]]+$SIGMA(VSUM(ev[i])/NCO(ev[i]))*1000.
     enddo
     
     sigma ns1=ns1+nskp
     sigma ns2=ns2+nskpi  
     
     sigma ens1=ens1+nskp*nskp
     sigma ens2=ens2+nskpi*nskpi
     
   enddo
   
   msig1=[msig1]/[nloop]
   msig2=[msig2]/[nloop]
   msig3=[msig3]/[nloop]
   
   st1=([msig1]+[msig2])/2.
   st2=([msig2]+[msig3])/2.
   
   message 'Average sigma k-p and k-pi: ' [st1] [st2]
   
   sigma evm=vmom/100.
   
   sigma ns1=ns1/[nloop]
   sigma ns2=ns2/[nloop]
   
   sigma ens1=SQRT((ens1-ns1*ns1*[nloop])/([nloop]-1.))
   sigma ens2=SQRT((ens2-ns2*ns2*[nloop])/([nloop]-1.))
   
   close 10
   fort/file 10 [outps]
   meta 10 -113
   
   title 'kaon separation'

   zone 1 2
   
* refractive indeces
* nC6F14 = 1.284
* nC5F12 = 1.259
   n=1.259
   
   message '*******'
   
*  pi mass
   m3=0.13957
* k mass
   m2=0.493677
* p mass
   m1=0.938
   
   do i=1,3
     j=40+[i]
     m=[m[i]]
     a[i]=acos(sqrt(1.+[m]*[m]/x/x)/[n])  
   enddo
   
   set txfp -60
   set chhe .5
   
   set pmci 1
   set mtyp 7
   set hcol 1

   null 2 5 0 50

*   graph [nbin] vmom ns1 paw
*   graph [nbin] vmom ns1 
   g/h/errors vmom ns1 evm ens1 [nbin] 20
   
   set fwid 10
   set fcol 2
   
* fun1 20 ([a2]-[a1])/([st1]/1000.) 1000 2. 10. s
   
   v/crea par(1) R 
*   v/fit vmom ns1 ens1 fnsig.f(x,[m2],[m1],[n])
   
   atitle ' ' 'k-p Separation (n[s])'
*   null 0 10 0 10 sab
   st1m=[st1]
   titolo=[s]?k-p! = [st1m] mr
* itx 5 7 [titolo]
   
   set pmci 1
   set mtyp 7
   set hcol 1

   null 2 5 0 20

*   graph [nbin] vmom ns2 paw
*   graph [nbin] vmom ns2 
   g/h/errors vmom ns2 evm ens2 [nbin] 22 
   
   set fcol 4
   
* fun1 20 ([a3]-[a2])/([st2]/1000.) 1000 2. 10. s
   
   atitle 'Momentum (GeV/c)' 'k-[p] Separation (n[s])'
*   null 0 10 0 10 sab
   st2m=[st2]
   titolo=[s]?k-[p]! = [st2m] mr
* itx 5 7 [titolo]
   set ltyp 3
   set plci 2
   line 2 4 5 4
   set plci 4
   line 2 3 5 3
   set plci 1
   set ltyp 1
   
   close 10
   
   close 1
   
   stringa=[wf] [surrad] [surdet] [msig1] [msig2] [msig3] [st1] [st2] $SIGMA(vnne(1)) $SIGMA(vnne(2)) $SIGMA(vnne(3))
   
   fmessage $QUOTE([stringa]) [outf]
   
   v/write vmom,ns1,ens1,ns2,ens2,v1,ev1,v2,ev2,v3,ev3,nv1,nv2,nv3 testout.txt '11f9.4,3f6.0' 

   message "############## Finito lrclas"
   
 return
 
* *******************************************************
* 
* 
 macro lim 1=1.4 2=1.7

   h/del *
   opt grid
   title 'Cherenkov and Limit Angles - [p], K, p at 2 GeV/c'
   zone
   null [1] [2] 0.4 1.
   fun1 110 asin(1./x) 100 [1] [2] sc
   for b in 0.9976 0.9709 0.9053
     fun1 120 acos(1./x/[b]) 100 [1] [2] sc
   endfor
   atitle 'n' '[q] (rad)'
 return
 
 
* *********************************************************
*
* 
 
 macro energy wf=1.0

   zone 1 1
   n/plot 10.1239.84/egamma/1e9 (z0g<[wf])&&(ipar==8)&&(1./(float(ngamma)+0.00001))
   
   atitle '[l] (nm)' 'E?[g]! distribution, electron (a.u.)'
   
 return
 
* *********************************************************

 macro qe 1=0

   title 'Photon Detector QE (from qe.dat)'
   v/del *
   
   v/read l,qe qe.dat
   v/read ll,qel qe-lambda.dat
   
   n=$VDIM(l,1)
   nn=$VDIM(ll,1)
   
   if ([1].eq.0) then
     sigma l=1239.84/l/1e9
   else
     sigma ll=1239.84/ll
   endif
   
   zone 1 2
   set mtyp 20
   set mscf .8
   graph [n] l qe pwc
   atitle ' ' 'QE in MC'
   graph [nn] ll qel pwc
   if ([1].eq.0) then
     atitle '[l] (nm)' 'qe boh'
   else
     atitle 'E (eV)' 'qe PC32-ALICE fig2.24'
   endif
   
 return

* **************************************************

macro optic 1=0

title 'Optical Characteristics (from optic.dat)'
v/del *
v/read l,fa,fe,fr,qa,qe,qr,ga,ge,gr,da,de,dr optic.dat
n=$VDIM(l,1)

* C5F12 refractive index

ea_u=21.364
fa_u=68.247
ga_u=([ea_u]**2-x**2)
fe_u=([fa_u]/[ga_u])
ff_u=sqrt((1.+2.*[fe_u])/(1.-[fe_u]))

sigma x = l*1e9
sigma frr = [ff_u]

v/print x
v/print frr

*

if ([1].eq.0) then
  sigma l=1239.84/l/1e9
endif

zone 2 4

tif='Freon C6F14'
tiq='Quarz'
tig='Gap (CH?4!)'
tid='Detector'

wif=1.
wiq=0.5
wig=10.
wid=.1
for j in f q g d
  sigma [j]b=exp(-[wi[j]]/([j]a+0.0001))
endfor

set ywin .2
set yval 9999.
set ylab .7
set mtyp 21
set mscf .6
for j in f q g d
 for i in b r
   if ([j].eq.'d') then
     set yval .1
   endif
   graph [n] l [j][i] pwl
   if ([j].eq.'d') then
     if ([i].eq.'b') then
       atitle 'Transparency - [l] (nm)' [ti[j]]
     else
	atitle 'Refractive Index - [l] (nm)'
     endif
   else
     if ([i].eq.'b') then
       atitle ' ' [ti[j]]
     endif
   endif
 endfor
endfor

return

*----------------------------------

macro trans

v/del ll,fft,qqt,ggt,e,ffa,qqa,gga
wiff=1.
wiqq=0.5
wigg=18.

set mtyp 6
opt grid

v/read ll,fft,qqt,ggt trans.dat
n=$VDIM(ll,1)

 title 'Absorbtion Length'
zone 1 3

sigma e=ll
*sigma e=1239.84/ll

for i in ff qq gg
  sigma [i]a=[wi[i]]/(log(1.)-LOG([i]t))
*  graph [n] ll [i]t pw
  graph [n] e [i]a pw
  atitle 'E (eV)' [i]
endfor

return

************************************************************

macro optra

exec proxie#optic
exec proxie#trans

n=$VDIM(l,1)
nn=$VDIM(ll,1)

zone 2 3

for i in f q g
  graph [n] l [i]a pw
  graph [nn] ll [i][i]a pw
endfor

return

************************************************************
* momentum dependence of angle and number of pe
*
macro vsmo wf=1.0 ef=1.0

zone 1 3

 opt utit
 set mtyp 6
 n/plot 10.prosel.f(1.,[wf],[ef])%sqrt(px**2+py**2+pz**2) prosel.f(1.,[wf],[ef])>0.

 atitle 'Momentum (GeV/c)' 'Cher Emission Angle (deg)'

 set mtyp 1
 n/plot 10.nfreon%sqrt(px**2+py**2+pz**2) 
 set mscf .2
 set mtyp 4
 n/plot 10.nfreon%sqrt(px**2+py**2+pz**2) ipar=8 option=s
 set mtyp 5
 n/plot 10.nfreon%sqrt(px**2+py**2+pz**2) ipar=11 option=s
 set mtyp 6
 n/plot 10.nfreon%sqrt(px**2+py**2+pz**2) ipar=14 option=s

 atitle 'Momentum (GeV/c)' 'Number of P.E.'

 zone 3 3 7 s 
 opt stat
 n/plot 10.nfreon ipar=8
 atitle '[p]' 'NPE'
 n/plot 10.nfreon ipar=11
 atitle 'K'
 n/plot 10.nfreon ipar=14
 atitle 'p'
 opt nsta

 opt htit
return

*******************************************************
* Compute the number of photons and photoelectrons
* Called by proxie#s
 
 macro n 1=8 ef=1.0 outf=' '

   par8='PI+'
   par11=_K+
   par14=__P

   m/gl/import nupe 0.
   m/gl/import nupad 0.

   zone 2 2

   opt stat
 
   n/plot 10.nfreon*[ef] ipar=[1]
   atitle 'N?PE!' 
   x1=$HINFO(1000000,'MEAN')
   x2=$HINFO(1000000,'RMS')
   x3=$HINFO(1000000,'ENTRIES')
   n/plot 10.nfreon ipar==[1].and.nfreon==0
   x0=$HINFO(1000000,'ENTRIES')
   
   n/plot 10.npad*[ef] ipar=[1]
   atitle 'N?PAD!'
   y1=$HINFO(1000000,'MEAN')
   y2=$HINFO(1000000,'RMS')
   y3=$HINFO(1000000,'ENTRIES')
   n/plot 10.npad ipar==[1].and.nfreon==0
   y0=$HINFO(1000000,'ENTRIES')
   
   a=Particle [par[1]]: [1] [x0] [x3] [x1] [x2] [y0] [y3] [y1] [y2]
   
   nupe=[x1]
   nupad=[y1]
   
   fmessage $QUOTE([a]) [outf]

 return

* **************************************************

 macro resume 1=1.

   opt utit

   titolo=Pa. Re.: Freon and gap widths scans, [e]=[1]
   title [titolo]
   h/del 20
   n/crea 20 'Resume' 8 ! ! wf wg tcro temi ttot npe npad effi
   n/read 20 resume.txt
   
   set ywin 1.2
   set ylab .6
   
   zone 2 2
   
   null 3 22 .5 7.
   set mtyp 20
   n/plot 20.ttot%wf wg=100.&&effi=[1] option=s
   set mtyp 28
   n/plot 20.temi%wf wg=100.&&effi=[1] option=s
   set mtyp 25
   n/plot 20.tcro%wf wg=100.&&effi=[1] option=s
   key 13. 3.5 20 '[s]?[q]!^c+e+p!'
   key 13. 3.0 28 '[s]?[q]!^c+e!'
   key 13. 2.5 25 '[s]?[q]!^c!'
   atitle 'w?freon! (w?g!=100) (mm)' '[s]?[q]! (mrad)'
   
   null 80. 240. .5 7.
   set mtyp 20
   n/plot 20.ttot%wg wf=14.&&effi=[1] option=s
   set mtyp 28
   n/plot 20.temi%wg wf=14.&&effi=[1] option=s
   set mtyp 25
   n/plot 20.tcro%wg wf=14.&&effi=[1] option=s
   key 100. 3.5 20 '[s]?[q]!^c+e+p!'
   key 100. 3.0 28 '[s]?[q]!^c+e!'
   key 100. 2.5 25 '[s]?[q]!^c!'
   atitle 'w?gap! (w?f!=14) (mm)' '[s]?[q]! (mrad)'
   
   null 3 22 5.*[1] 40.*[1]
   set mtyp 21
   n/plot 20.npe*effi%wf wg=100.&&effi=[1] option=s
   set mtyp 26
   n/plot 20.npad*effi%wf wg=100.&&effi=[1] option=s
   key 15. 15.*[1] 21 'N?P.E.!'
   key 15. 10.*[1] 26 'N?PAD!'
   atitle 'w?freon! (w?g!=100) (mm)' 'N'
   
   null 80. 240. 10.*[1] 28.*[1]
   set mtyp 21
   n/plot 20.npe*effi%wg wf=14.&&effi=[1] option=s
   set mtyp 26
   n/plot 20.npad*effi%wg wf=14.&&effi=[1] option=s
   key 100. 16.*[1] 21 'N?P.E.!'
   key 100. 14.*[1] 26 'N?PAD!'
   atitle 'w?gap! (w?f!=14) (mm)' 'N'
   
 return
 
* **************************************************
* Draw ....
* 
* mip distribution
* fired pad distribution
*  open the rz file with: proxie nomefile
* phthr = number of pad threshold (for efficiency estimation)
* rot: if 1 the pad are rotated respect to the hypernuclear experiment
* 
* r1f = smaller radius of radiator OR radiator definition file
* r1d = smaller radius of detector OR detector definition file
 
 macro ophase 1=proxie phthr=5 r1f=49. r2f=265. xof=-222.5 anf=27.1 r1d=49. r2d=357. xod=-269. a1d=-60. a2d=60. ofile='paw_ophase.eps'
   
   opt ngri
   close 1
   
   ifile=[1].rz
   set lwid 1 
   IF ($FEXIST([ifile]) .EQ. 0) THEN
     MESSAGE [1] file does not exist
     GOTO fine_ophase
   ENDIF
   
   exec proxie [1].rz
   
   IF ($FEXIST([r1f]) .EQ. 1) THEN 
     exec proxie#loardg [r1f]
     
     r1f=$SIGMA(vrdg(1))
     r2f=$SIGMA(vrdg(2))
     xof=$SIGMA(vrdg(3))
     anf=$SIGMA(vrdg(4))
     
   ENDIF
    
   IF ($FEXIST([r1d]) .EQ. 1) THEN 
     exec proxie#loaddg [r1d]
     
     r1d=$SIGMA(vddg(1))
     r2d=$SIGMA(vddg(2))
     xod=$SIGMA(vddg(3))
     a1d=$SIGMA(vddg(4))
     a2d=$SIGMA(vddg(5))
     
   ENDIF
   
* compute acceptance efficiency
   
    zone 1 2
    n/plot 10.y0%x0 (ipar.eq.11).and.(npad.le.[phthr])
    n0k=$HINFO(1000000,'ENTRIES')
    n/plot 10.y0%x0 (ipar.eq.11)
    nnk=$HINFO(1000000,'ENTRIES')
 
    eff=([nnk]-[n0k])/[nnk]
    message 'Charge Particle with <= ' [phthr] ' photons:' [n0k]
    message 'Charge Particle with any photon(s):' [nnk]
    message 'k Acceptance Efficiency: ' [n0k] [nnk] [eff]

* plot phase space
   
    title [1]
    opt nsta
   
    close 10
    fort/file 10 [ofile]
    meta 10 -113
   
    zone 1 1

   null -400 400 -400 400 

    set mtyp 1
    set pmci 7
* n/plot 10.yg%xg (yg.le.100)&&(xg.le.100)&&(yg.ge--100)&&(xg.ge.-200) option=cont
    n/plot 10.yg%xg prosel.f(6.,0.,0.) option=s
   
* set pmci 5
    set mtyp 20
    set mscf .5
	set pmci 2
    n/plot 10.ygp%xgp prosel.f(6.,0.,0.) option=s
    atitle 'x (cm)' 'y (cm)'
    set pmci 1
    n/plot 10.y0%x0 prosel.f(6.,0.,0.) option=s
    
* radiator
    exec proxie#drawarc ! col=4 r1=[r1f] r2=[r2f] x0=[xof] y0=0. an1=-[anf] an2=[anf]
* detector
    exec proxie#drawarc ! col=2 r1=[r1d] r2=[r2d] x0=[xod] y0=0. an1=[a1d] an2=[a2d]
    
    null 0 10 0 10 sab
    
* itx 1 9  'MIP (red) and PE/HITS Spatial Distribution'
    testo=k Acc. Efficiency = $FORMAT([eff],F5.3) (pad.gt.[phthr])
    itx .5 9.5 [testo]
    close 10
    
    message '---------------------'
    message '-- Output on: ' [ofile]
    message '---------------------'
fine_ophase:
 
 
 return
 
* **************************************************
* 
* Plot the particle (trajectories) phase space 
* at RICH entrance window
 
 macro phase 1=0 2='none'
      
   opt utit
   title 'Input Particle Phase Space'
   set mtyp 6 
   set ywin 1.5
   
   if ([1].eq.0) then
     zone 3 2
     n/plot 10.thetapar%phipar phipar<-2.
     atitle ' ' 'theta (rad)'
   else
     zone 1 2
   endif
   null -0.05 0.05 -0.01 0.15
   n/plot 10.thetapar%phipar -2.<phipar<2.  option=s
   if ([1].eq.0) then 
     n/plot 10.thetapar%phipar 2.<phipar
     atitle 'phi (rad)'
   else
     atitle 'phi (rad)' 'theta (rad)'
   endif
   
   zone 2 2 3 s
   null -110 110 -16 16
   n/plot 10.y0%x0  option=s
   atitle 'x?0! (cm)' 'y?0! (cm)'
   n/plot 10.sqrt(px**2+py**2+pz**2)
   atitle 'Momentum distribution (GeV/c)'
 return
 
* ***************************
* 
 
 macro drawarc opt='s' col=2 r1=49. r2=265. x0=-222.5 y0=0. an1=-27.1 an2=27.1
   
   set plci [col]
   set lwid 5
   set bord 1
   g/pri/arc [x0] [y0] [r1] [r2] [an1] [an2]
   
 return
 
* ***********************************************************
* 
* Draw the radiator vessel
 
 macro drawrad opt='s' col=3 wid=5
   
* unit in cm
   
   rad_sx=91.05
   rad_sy=15.5
   
   xmax=[rad_sx]*1.1
   ymax=[rad_sy]*1.1
   
   if ([opt] .ne. 's') then
     null -[xmax] [xmax] -[ymax] [ymax]
   endif
   
   set plci [col]
   set lwid [wid]
   box -[rad_sx] [rad_sx] -[rad_sy] [rad_sy]
   set lwid 1
   
   
 return
 
 
 * ***********************************************************
* 
* Draw the pad planes
* if pad_gy = 0 draw 3x1 planes, 3x2 otherwise
* 
* if opt = s overlay on top of existing picture
* 
* nx = pad along x
* ny = pad along y
* pad_sx, pad_sy = semisize of the pad plane
* pad_gx, pad_gy = distance bettwen pad planes
* col, wid : color and width fo the lines
* s: if set overwrite on top of existing figure
 
 macro drawpad pad_gx=2.0 pad_gy=0 pad_sx=32. pad_sy=20.16 nx=3 ny=1 col=4 wid=5 opt='s'
   
   set lwid [wid]
   set pmci [col]
   
   sx=([pad_sx]*2.+[pad_gx])
   sy=([pad_sy]*2.+[pad_gy])
   
   dx=[pad_gx]
   dy=[pad_gy]
   
   xmax=[sx]*[nx]/2.
   ymax=[sy]*[ny]/2.
	
   if ([opt] .ne. 's') then
     xxmax=1.1*[xmax]
     yymax=1.1*[ymax]
     null -[xxmax] [xxmax] -[yymax] [yymax]
   endif

   do j=1,[ny]
     my = -[ymax]+[sy]*([j]-1.)
     py0=[my]+[dy]/2.
     py1=[my]+[sy]-[dy]/2.

     do i=1,[nx]
       mx = -[xmax]+[sx]*([i]-1.)
       
       px0=[mx]+[dx]/2.
       px1=[mx]+[sx]-[dx]/2.
       

       box [px0] [px1] [py0] [py1]
     
     enddo
   enddo
   
   set lwid 1
   
 return
 
 
* ***********************************************************
* 
* Draw the pad planes
* if pad_gy = 0 draw 3x1 planes, 3x2 otherwise
* 
* if opt = s overlay on top of existing picture
* 
* nx = pad along x
* ny = pad along y
* rot = if 1, pad are rotated 90 degree (sx<sy)
 
 macro drawpadold pad_gy=0 opt='s' col=4 wid=5 
   
   pad_sx = 32.
   pad_sy = 20.16
   
   IF ([rot] .EQ. 1) THEN
     dummy=[pad_sx]
     pad_sx = [pad_sy]
     pad_sy = [pad_sx]
   ENDIF
   
   pad_gx = 2.
   
   set lwid [wid]
   set pmci [col]
 
   pad_x = 3.*[pad_sx]+[pad_gx]
   pad_x1 = [pad_sx]+[pad_gx]
   pad_x2 = [pad_sx]
   
   xmax=1.05*[pad_x]
   
   if ([pad_gy] .eq. 0.) then
     
     pad_y = [pad_sy]
     ymax = 1.1*[pad_y]
     
     if ([opt] .ne. 's') then
       null -[xmax] [xmax] -[ymax] [ymax]
     endif
     
     set plci [col]
     
     box -[pad_x] -[pad_x1] -[pad_y] [pad_y]
     box -[pad_x2] [pad_x2] -[pad_y] [pad_y]
     box [pad_x] [pad_x1] -[pad_y] [pad_y]
     
   else
   
     pad_y = 2.*[pad_sy]+[pad_gy]/2.
     pad_y1 = [pad_gy]/2.
     ymax = 1.05*[pad_y]
     
     if ([opt] .ne. 's') then
       null -[xmax] [xmax] -[ymax] [ymax]
     endif
     
     set plci [col]
   
     box -[pad_x] -[pad_x1] -[pad_y] -[pad_y1]
     box -[pad_x] -[pad_x1] [pad_y] [pad_y1]
     box -[pad_x2] [pad_x2] -[pad_y] -[pad_y1]
     box -[pad_x2] [pad_x2] [pad_y] [pad_y1]
     box [pad_x] [pad_x1] -[pad_y] -[pad_y1]
     box [pad_x] [pad_x1] [pad_y] [pad_y1]
   
   endif
   
   set lwid 1
   
 return
 
 
***************************************************
* 
* Plot a single event
* pad_gy: 0 -> 3 plane configuration, >0 -> six plane configuration
 
 macro single 1=1 wf=1.5 pad_gy=0
   
   titolo=single event [1]
   title [titolo]
* zone 1 4
 zone 1 1
   opt nsta
   opt utit
   set ywin 0.6

   set pmci 3
 
   exec proxie#drawpadold [pad_gy] opt='nnn' col=1
 
   set mtyp 23
   set pmci 3
   n/plot 10.yg%xg ivent=[1] option=s
   
   set mtyp 21
   set pmci 4
   n/plot 10.ygp%xgp ivent=[1] option=s
   
   set mtyp 26
   set pmci 2
   n/plot 10.y0%x0 ivent=[1] option=s

   message 'return to continue'
   read dummy

   opt sta
   set hcol 1
   n/plot 10.theta0g/180.*3.1415 ivent=[1]

   set hcol 2
   n/plot 10.prosel.f(1.0,[wf],1.0)/180.*3.1415 ivent=[1] option=s

   set hcol 1
   n/plot 10.etha ivent=[1]

   set hcol 2
   n/plot 10.prosel.f(2.0,[wf],1.0) ivent=[1] option=s

   set hcol 1
   n/plot 10.rthpad ivent=[1]

   set hcol 2
   n/plot 10.prosel.f(3.0,[wf],1.0) ivent=[1] option=s
   
   set hcol 1
   
 return

* **************************************************
* Phase Space from Optics (LeRose) File *.mud
* 
* dz is the distance (in meter) along the TRANSPORT System z
* 
* NOTE: mud data of the hadron arm refer to the rich entrance window at about 2.6 m 
* from the TRANSPORT system origin (on the center of the bottom VDC)
* NOTE: the electron mud parameters refer to the origin of the TRANSPORT system
* 

 macro phace 1='rich.mud' dz=0.0
  
   close 1
   
   close 10
   fort/file 10 paw_phace.eps
   meta 10 -113
   
   opt utit

   set xlab 1.2
   set ylab .7
   set ywin 1.2
   
   set vsiz .3
   set asiz .4
   
   n/crea 123 'PhaseSpace' 10 ! ! x0 t0 y0 p0 xf tf yf pf l dp
   
   n/read 123 [1]
   
   message 'Data read from file: ' [1]
   message 'Data now in ntuple: ' 123

   zone 3 3
   
   ddz=2.6+[dz]
   titolo=Trajectory Phase Space at z=[ddz] m
   title [titolo]
   
   n/plot 123.(xf+tan(tf)*[dz])
   atitle 'X (m)'

   n/plot 123.(yf+tan(pf)*[dz])
   atitle 'Y (m)'

   n/plot 123.dp
   atitle 'dp/p'
   
   n/plot 123.pf
   atitle '[f]'
   
   n/plot 123.tf
   atitle '[q]'
   
   n/plot 123.(yf+tan(pf)*[dz])%(xf+tan(tf)*[dz])
   atitle 'X (m)' 'Y (m)'
   
   n/plot 123.pf%tf
   atitle '[q]' '[f]'

   n/plot 123.tf%(xf+tan(tf)*[dz])
   atitle 'X (m)' '[q]'

   n/plot 123.pf%(xf+tan(tf)*[dz])
   atitle 'X (m)' '[f]'
   
   close 10
   
 return

* **********************************************************************
*
* read summary text as: thick err1_k err1_p err2_k err2_p err3_k err3_p
*
* err1 = emission
* err2 = + pint pad
* err3 = + finite pad
*
 macro rsu 1='summary_lerose_65.txt'

title 'RICH Gap (LeRose optics input, [e]=0.65)'
set yval 999.
set ywin .2
set xlab 1.2

v/read t,a1,b1,a2,b2,a3,b3 [1]

xa1=1.75
xa2=1.
xa3=1.
xb1=2.3
xb2=1.
xb3=1.

size 20 20
zone 1 2

n=$VDIM(t,1)

sigma a1=a1*[xa1]
sigma a2=a2*[xa2]
sigma b1=b1*[xb1]
sigma b2=b2*[xb2]

message [n]

m=3*[n]

v/crea xx([m])
v/crea ya([m])
v/crea yb([m])

j=1
do i=1,3 
v/copy t xx([j]:)
v/copy a[i] ya([j]:)
v/copy b[i] yb([j]:)
j=[j]+[n]

enddo

v/crea par(3)
v/crea minx(6)

k=1
tm = $SIGMA(t([n]))

for i in a b

  set mtyp 1
  graph [m] xx y[i] paw
  jj=20

  do j=1,3
    set mtyp [jj]
    graph [n] t [i][j] p
    sigma y=[i][j]/100.
    v/fit t [i][j] y p2 S 3 par

    p1=$SIGMA(par(1))
    p2=$SIGMA(par(2))
    p3=$SIGMA(par(3))
    ym=$SIGMA([p1]+[p2]*[tm]+[p3]*[tm]*[tm])

    mx=$SIGMA(-par(2)/2./par(3))

    tmx=[tm]+.75

*    testo="X#[x[i][j]]
    testo=$FORMAT([mx],F5.2)
    itx [tmx] [ym] [testo]

    message [tm] [ym] [x[i][j]]
    v/input minx([k]:[k]) [mx]
    k=[k]+1
    jj=[jj]+2
  enddo

  if ([i] .EQ. 'a') then 
     key 13 7.4 20 'Emission'
     key 13 7.8 22 'Point Pad'
     key 13 8.2 24 'Finite Pad'

    set yval .15
    atitle '' '[s] Kaon (mr)'
  endif

endfor

atitle 'Gap (cm)' '[s] Pion (mr)'

v/pri minx

return

* -----------------------------------------------------------
*
* read new summary text file
*
* out_c5f12_gap_2gev.txt

macro rsun 1='out_c6f14_gap.txt'
  
  set *fon -60
  set chhe .5
  set gsiz .45
  set csiz .4
  set asiz .4
  set vsiz .4
  set ygti 1.
  option grid
  option nbox
  
  set yval 999.
  set ywin .2
  set xlab 1.2
 
 
  zone 1 2

  v/del *

* title 'C6F14 - LeRose optics - 2 GeV/c - 12.5 PADs/10 cm'
  title 'C5F12 - LeRose optics - 2 GeV/c'
  title 'C6F14 - LeRose optics - 3 GeV/c'
  title 'C6F14 - 2001 phase space - 2 GeV/c'
  title 'C5F12 - no septum - 2 GeV/c'
  title 'C6F14 - no septum - 3 GeV/c'
  title 'C6F14 - septum - 3 GeV/c'
  title 'C6F14 - septum - 2 GeV/c'
  title 'C?6!F?14!, 2.28-2.52 GeV/c, 2.0/6 planes, 8x8.4 mm^2!, z=2.6 m'
  title 'C?6!F?14!, 2.28-2.52 GeV/c, 3.5/6 planes, 8x8.4 mm^2!, z=2.6 m'
  title 'C?6!F?14!, 2.28-2.52 GeV/c, 3 planes, 8x8.4 mm^2!, z=2.3 m'
  title 'C?6!F?14!, 2.28-2.52 GeV/c, 3.5/6 planes, 8x8.4 mm^2!, z=2.3 m'
  
  title 'C?6!F?14!, 2.28-2.52 GeV/c, 3.5/6 planes, 8x8.4 mm^2!, z=1.8 m'
 
  title [1]
 
  gap_0 = 8.
  gap_1 = 36.
 
  v/read t,ns1,ns2,ns3,sk1,sk2,sk3,spi1,spi2,spi3,pep,pap,pek,pak,pepi,papi [1]

  v/crea par(3)

  n=$VDIM(t,1)
  
  sigma xt=t*2.

* null 5. 21 4 20

  ma1=$SIGMA(VMAX(ns1))
  ma2=$SIGMA(VMAX(ns2))
  ma3=$SIGMA(VMAX(ns3))
  mi1=$SIGMA(VMIN(ns1))
  mi2=$SIGMA(VMIN(ns2))
  mi3=$SIGMA(VMIN(ns3))

  f1=[ma1]/[ma3]*.7
  f2=[ma2]/[ma3]*.9
 
  f1=1.5
  f2=1.
 
  ma1=[ma1]/[f1]
  ma2=[ma2]/[f2]
  mi1=[mi1]/[f1]
  mi2=[mi2]/[f2]
  
  message [ma1] [ma2] [ma3] [mi1] [mi2] [mi3]
  
  ma1=10.
  ma2=[ma1]
  ma3=[ma1]
  mi1=4.
  mi2=[mi1]
  mi3=[mi1]
  
  v/crea mai(6) R [ma1] [ma2] [ma3] [mi1] [mi2] [mi3]
  
  sigma ns1=ns1/[f1]
  sigma ns2=ns2/[f2]
  
  sigma dy=ns1/20.
  sigma dx=t/t*.001
  
  close 10
  
  fort/file 10 paw_out_rsun.eps
  meta 10 -113
  
  y2=$SIGMA(VMAX(mai))*1.1
  y1=$SIGMA(VMIN(mai))*.9
  
  null [gap_0] [gap_1] [y1] [y2]
  
  set mtyp 20
  set pmci 5
  set plci 5
* graph [n] xt ns1 p
  sigma dy=ns1*.025
  graph/h/errors xt ns1 dx dy [n] 20 .2 
* v/fit xt ns1 dy p2 SW 3 par	
  set mtyp 25
  set pmci 4
  set plci 4
* graph [n] xt ns2 p
  sigma dy=ns2/20.
  graph/h/errors xt ns2 dx dy [n] 25 .2 
* v/fit xt ns2 dy p2 WS 3 par
  set mtyp 22
  set pmci 6
  set plci 6
* graph [n] xt ns3 p
  sigma dy=ns3*.025
  graph/h/errors xt ns3 dx dy [n] 22 .2 
* v/fit xt ns3 dy p2 WS 3 par
  
  atitle '' 'n?[s]!(k,[p])'
  
  null 0 10 0 10 sab
  
  key 6.5 8.2 22 'Finite Pad'
  set pmci 4
  yt=("X#$FORMAT([f2],F3.1)) Point Pad
  key 6.5 8.8 25 [yt]
  set pmci 5
  yt=("X#$FORMAT([f1],F3.1)) Emission Only
  key 6.5 9.5 20 [yt]
  
  set yval .2
  
* zone 1 1
  
  null [gap_0] [gap_1] 8.0 18.
  
  npappa=.8
  
  sigma pap=pap/[npappa]
  
  set mtyp 20
  set pmci 1
  graph [n] xt papi p
  set mtyp 22
  set pmci 2
  graph [n] xt pak p
  
  set mtyp 24
  set pmci 3
  graph [n] xt pap p
  
  null 0 10 0 10 sab
  testo=("X#[npappa]) p
  key 6.5 8.2 24 [testo]
  set pmci 2
  key 6.5 8.8 22 'K+'
  set pmci 1
  key 6.5 9.5 20 '[p]+'
  
  atitle 'Gap (cm)' 'Mean number of firing PADs'
  
  close 10

return

****
*
*
macro y

sigma t=2.*t

xa1=1.75
xa2=1.
xa3=1.
xb1=2.3
xb2=1.
xb3=1.

size 20 20
zone 1 2

n=$VDIM(t,1)

sigma a1=a1*[xa1]
sigma a2=a2*[xa2]
sigma b1=b1*[xb1]
sigma b2=b2*[xb2]

message [n]

m=3*[n]

v/crea xx([m])
v/crea ya([m])
v/crea yb([m])

j=1
do i=1,3 
v/copy t xx([j]:)
v/copy a[i] ya([j]:)
v/copy b[i] yb([j]:)
j=[j]+[n]

enddo

v/crea par(3)
v/crea minx(6)

k=1
tm = $SIGMA(t([n]))

for i in a b

  set mtyp 1
  graph [m] xx y[i] paw
  jj=20

  do j=1,3
    set mtyp [jj]
    graph [n] t [i][j] p
    sigma y=[i][j]/100.
    v/fit t [i][j] y p2 S 3 par

    p1=$SIGMA(par(1))
    p2=$SIGMA(par(2))
    p3=$SIGMA(par(3))
    ym=$SIGMA([p1]+[p2]*[tm]+[p3]*[tm]*[tm])

    mx=$SIGMA(-par(2)/2./par(3))

    tmx=[tm]+.75

*    testo="X#[x[i][j]]
    testo=$FORMAT([mx],F5.2)
    itx [tmx] [ym] [testo]

    message [tm] [ym] [x[i][j]]
    v/input minx([k]:[k]) [mx]
    k=[k]+1
    jj=[jj]+2
  enddo

  if ([i] .EQ. 'a') then 
     key 13 7.4 20 'Emission'
     key 13 7.8 22 'Point Pad'
     key 13 8.2 24 'Finite Pad'

    set yval .15
    atitle '' '[s] Kaon (mr)'
  endif

endfor

atitle 'Gap (cm)' '[s] Proton (mr)'

v/pri minx

return

******************************************
* C6F14 and C5F12 Refractive Indices
*

macro nfreon

zone 1 1
option utit

*C5F12

ea=21.364
fa=68.247
ga=([ea]**2-x**2)
fe=([fa]/[ga])

ff12=sqrt((1.+2.*[fe])/(1.-[fe]))

fun1 10 [ff12] 100 5. 11.

*C6F14

a=1.1777
b=0.0172

fn=[a]+[b]*x

fun1 11 [fn] 100 5. 11. s

return

endkumac