macro make_figure pt1=0.0 pt2=0.3 pt3=0.6 pt4=0.9
*
* this is a specialized, one-time-only macro to make a plot showing
*  1) the HBT radii oscillate
*  2) it is a second-order oscillation
*
*
vec/crea ptcuts(4) r [pt1] [pt2] [pt3] [pt4]
*
op ndat
op utit
set lwid 8
set fwid 8
set hwid 8
set csiz 0.5
set vsiz 0.6
set xval 0.45
set ylab 0.9
set xlab 2.1
set xwin 0.001
set ywin 0.001
set asiz 0.55
set bcol 1
set chhe 0.55
set txal 20
set xmgr 2.7
*
size 20 35
2dhisto 1 'Ro2'  40 0.0 6.4 40 5.1 39.9
2dhisto 2 'Rs2'  40 0.0 6.4 40 5.1 39.9
2dhisto 3 'Ros2' 40 0.0 6.4 40 -8 8
2dhisto 4 'Rl2'  40 0.01 6.4 40 5 130
zone 2 2
exec make_figure#plot_histo 1
set xlab 1.5
*atit ' ' 'R^2?O! (fm^2!)' ! 2
set TANG 90
itx -1.4 30.0 'R^2?O! (fm^2!)'
*
exec make_figure#plot_histo 2
set xlab 1.9
*atit ' ' 'R^2?S! (fm^2!)' ! 2 R
set TANG 270
itx 8.1 30.0 'R^2?S! (fm^2!)'
*
exec make_figure#plot_histo 3
set xlab 1.5
atit '[f]?p! (rad) '
set TANG 90
itx -1.4 4.0 'R^2?OS! (fm^2!)'
*
exec make_figure#plot_histo 4
set xlab 1.9
atit '[f]?p! (rad) '
set TANG 270
itx 8.1 100.0 'R^2?l! (fm^2!)'
*

exec make_figure#read_data
*
first = 1

ykey = 90
do ipt=1,4
  set basl 0.01
  case [ipt] in
  (1) lt = 1
  (2) lt = 12
  (3) set basl 0.02
      lt = 12
  (4) lt = 13
  endcase
  exec make_figure#sort_data [ipt]
  do ipar=1,4
*    trace on
*    mess ipar = [ipar]
*    zone 2 2 [ipar] S
* paw is a fucking piece of shit!!
    case [ipar] in
    (1) graphics/viewing/zone 2 2 1 s
	exec make_figure#plot_histo 1
    igset ltyp [lt]
    set ltyp [lt]
        pline $VDIM(phival) phival Ro2_sel
    (2) graphics/viewing/zone 2 2 2 s
	exec make_figure#plot_histo 2
    igset ltyp [lt]
    set ltyp [lt]
        pline $VDIM(phival) phival Rs2_sel
    (3) graphics/viewing/zone 2 2 3 s
	exec make_figure#plot_histo 3
    igset ltyp [lt]
    set ltyp [lt]
        pline $VDIM(phival) phival Ros2_sel
    (4) graphics/viewing/zone 2 2 4 s
	exec make_figure#plot_histo 4
    igset ltyp [lt]
    set ltyp [lt]
        pline $VDIM(phival) phival Rl2_sel
    endcase
*    trace off
  enddo

 igset ltyp [lt]
 set ltyp [lt]
 ptvalue = $SIGMA(ptcuts([ipt]))
 key 1.3 [ykey] [lt] 'p?T! = '//[ptvalue]//' GeV/c' 1.2 L
 ykey = [ykey]-10
enddo
return

macro plot_histo ih
case [ih] in
(1) set yval 200.
    set xval 0.3
    set ndvy 10510.00
(2) set yval 200.
    set xval -7.7
    set ndvy 10510.01
(3) set yval 0.2
    set xval 0.3
    set ndvy 10510.00
(4) set yval 0.2
    set xval -7.7
    set ndvy 10510.01
endcase
hi/plo [ih]
tic
return



*
*-----------------------------------------------------------
macro read_data 
*
vec/rea pt,phi_p,N,x2,y2,xy,xt,yt,t2,z2,x,y,t,Ro2,Rs2,Ros2,Rl2 Calc1/pions-vs-phi.BW ! ! -/%/
return
*
*-----------------------------------------------------------


*-----------------------------------------------------------
macro sort_data [ipt]
ptcutval = $SIGMA(ptcuts([ipt]))
*
* count how many elements I need...
ict = 0
do i=1,$VDIM(pt)
  if ($SIGMA(pt([i])).eq.[ptcutval]) then
    ict = [ict]+1
  endif
enddo

vec/del  Ro2_sel
vec/crea Ro2_sel([ict]) r
vec/del  Rs2_sel
vec/crea Rs2_sel([ict]) r
vec/del  Ros2_sel
vec/crea Ros2_sel([ict]) r
vec/del  Rl2_sel
vec/crea Rl2_sel([ict]) r
vec/del  phival
vec/crea phival([ict]) r
* now fill 'em
ict = 0
do i=1,$VDIM(pt)
  if ($SIGMA(pt([i])).eq.[ptcutval]) then
    ict = [ict]+1
    vec/in phival([ict])   $SIGMA(phi_p([i]))
    vec/in Ro2_sel([ict])  $SIGMA(Ro2([i]))
    vec/in Rs2_sel([ict])  $SIGMA(Rs2([i]))
    vec/in Ros2_sel([ict]) $SIGMA(Ros2([i]))
    vec/in Rl2_sel([ict])  $SIGMA(Rl2([i]))
  endif
enddo
return
*-----------------------------------------------------------