macro digi igo=1
* sh  igo=1 -> P0 ; igo=2 -> P0'
if [igo]=1 then
*P0
  p = 'p0_2_1_20000'
  pn = 'p0'
else
*P0'
  p = 'p0_3_1_20000'
  pn = 'p0p'
endif
hi/file 1 [p].rz 0; nt/print  400; call digi.f(0); 
close 1
* sh makes an output file of all the histograms
hi/file 1 [pn].rz ! n
hrout *
close 1
return
**************************************************
macro rhit igo=1
if [igo]=1 then
*P0
  p = 'p0_2_1_20000'
else
*P0'
  p = 'p0_3_1_20000'
endif
hi/file 1 [p].rz 0; nt/print 900; hrin 900; call rhit.f; close 1
return
**************************************************
macro plot igo=1 opt=depe ch=42 dd=405 scn=1
opt utit
opt stat
if [igo]=1 then
  www = 'p0/'
*  www = './'
  det = 'p0'
else 
  www = 'p0p/'
  det = 'p0p'
endif
fn = [www][opt].ps
if [opt]=depe then
  fn = [www]depenergy.ps
endif
if [opt]=sigma then
  if [scn]=2 then
    fn = [www]sigmas2n.ps
  endif
endif
zone
for/file 67 [fn]
meta 67 -114
next
case [opt] in
    (depe)   zone 1 3 ; 
             set stat [ch]100; opt liny
             title 'Deposited energy' ; 
             hi/pl [ch]001; 
             atitle 'Deposited energy (keV)' ' ' ; 
             hi/pl [ch]100; 
             atitle '[S] ADC' ''
             hi/pl [ch]702; 
             atitle 'Signal-to-Noise Ratio' ''
    (stn)    title 'Signal-to-Noise Ratio for P0 Test Chamber,Simulated vs. Real Data' ; 
             hi/pl [ch]702; 
             atitle 'Signal-to-Noise Ratio' ''
*
    (lr)     opt nsta; igset mtyp 1; 
	     title_g 'Right versus Left for lower corner of chamber';
             hi/pl [ch]200; 
             atitle 'Left to Center' 'Right to Center'
*
    (lrw)    opt nsta; igset mtyp 1; 
	     title_g 'Right versus Left for only wire group'; 
             hi/pl [ch]210 ; 
             atitle 'Left to Center' 'Right to Center'
*
    (dphi)   opt nsta; igset mtyp 1; 
	     title 'Distance from right upper corner versus [F]'
             hi/pl [ch]201; atitle '[F] (degrees)' 'Distance'
*
    (dphiw) opt nsta; igset mtyp 1; 
	    title 'Distance from right upper corner versus Phi for only wire group'
           hi/pl [ch]211; atitle '[F] (degrees)' 'Distance'
*
    (profdphi) opt nsta; igset mtyp 20; 
	       title 'Distance from right upper corner versus [F]'
          min [ch]202 0.8;  hi/pl [ch]202; atitle '[F] (degrees)' 'Distance'
*
    (profdphiw) opt nsta; igset mtyp 20; 
	        title 'Distance from right upper corner versus Phi for only wire group'
           min [ch]212 0.8;  hi/pl [ch]212; atitle '[F] (degrees)' 'Distance'
*
    (vlad) opt nsta; igset mtyp 1;    
       	   title 'Relative charge distribution'; 
           hi/pl [ch]221 ; atitle 'X (strip)' 'Relative charge'
*
    (vladw) opt nsta; igset mtyp 1;   
	    title 'Relative charge distribution for only wire group'
            hi/pl [ch]222 ; atitle 'X (strip)' 'Relative charge'
*
    (res) id=[ch]303;   opt stat; igset mtyp 20;    
	  title 'Resolution (with known wire)'; 
	  opt fit; hi/fit [id](-0.02:0.01) g  0;
          hi/pl [id](-0.05:0.05); atitle '[D]X (cm)' ''
*
    (ress) id=[ch]301; opt stat; igset mtyp 20;    
	   title 'Resolution'; 
	   opt fit; hi/fit [id](-0.05:0.05) g  0;
           hi/pl [id](-0.15:0.15); atitle '[D]X (cm)' ''
*
    (resp)  opt nsta; igset mtyp 20;    
    	    title 'Resolution versus position on strip'; 
	    opt nsta; zone 1 2
               for id in [ch]402 [ch]404 
                 min [id] -100.; max [id] 100; hi/pl [id]
                 null 0 1 0 1 sab
                 if (mod([id],10)=2) then
                   text 0.2 0.9 'All' 0.4 
*   Need to adjust microns back to cm when digitization is fixed
                   atitle ' ' '[D]X ([m]m)'
                 elseif (mod([id],10)=4) then
                   text 0.2 0.9 'for only wire group' 0.4 
                   atitle 'X (cm)' ' '
		 endif
               endfor
               igset mtyp 1
*
    (sigma)  opt nsta; igset mtyp 1; 
	     title 'Resolution (sigma) versus position on strip'; 
	     if [scn]=2 then 
	       title 'P0 Signal-to-Noise Ratio, Resolution vs. strip position'; 
	     endif
	     zone 1 2 
	     opt nsta
             nbins = 12
             vec/cre par(3); vec/cre err(3); vec/cre av([nbins]); 
	     vec/cre eav([nbins])
             vec/cre sigmf([nbins]); vec/cre esigmf([nbins])
             id =[ch][dd]
             id2 = [id]-2
             igset mtyp 1
             sliy [id] [nbins]
             hi/proj [id]
             opt stat
             opt fit
             nx  = $HINFO([id],'XBINS')
             xmn = $HINFO([id],'XMIN')
             xmx = $HINFO([id],'XMAX')
             dx  = ([xmx]-[xmn])/[nx];
             sigma xbin = array([nbins],[xmn]+[dx]/2#[xmx]-[dx]/2)
        do ib = 1, [nbins]
          hi/cop [id].sliy.[ib] 99  
*          hi/fit 99 g ! 0 par ? ? ? err
          hi/fit 99 g 0 0 par ? ? ? err
*             rms =   $HINFO(99,'RMS')
*             mean =  $HINFO(99,'MEAN')
*             xxmin = [mean]-2*[rms]
*             xxmax = [mean]+2*[rms]
*          hi/fit 99([xxmin]:[xxmax]) g ! 0 par ? ? ? err
*          wait
          hi/fit 99([xxmin]:[xxmax]) g 0 0 par ? ? ? err
          vect/copy par(2) av([ib])
          vect/copy err(2) eav([ib])
          vect/copy par(3) sigmf([ib])
          vect/copy err(3) esigmf([ib])
        enddo
        igset mtyp 1; opt nsta; opt ndat; 
	if [scn]=1 then
          hi/pl [id](,-2500.:2500.)  
          atitle 'relative position on strip' 'Residual ([m]m)'
	else if [scn]=2 then
	  hi/pl 42701
	  atitle 'Signal-to-Noise Ratio'
	endif
        igset mtyp 20
        sigma av = av; sigma eav = eav
        sigma sigmf = sigmf; sigma esigmf = esigmf
* sh  Plots average distribution of [id] in [nbins]
*        null  [xmn] [xmx] -100. 100.
*	hpl/err xbin av    ? eav    [nbins] 20 0.25
*        atitle 'relative position on strip' 'Average ([m]m)'
* sh  Plots sigma vs. delx
	if ([ch] = 10)  then 
	  ylim = 250.
	elseif ([ch] = 11)  then 
	  ylim = 250.
	elseif ([ch] > 11)  then 
	  if [igo]=1 then
	    ylim = 500. 
	  else 
	    ylim = 100.
	  endif
	endif
	null  [xmn] [xmx] 0. [ylim]  
        hpl/err xbin sigmf ? esigmf [nbins] 20 0.25
	key -0.45 0.9*[ylim] 20 'Simulation'
        atitle 'relative position on strip' '[s] ([m]m)'
* sh  This plots the experimental data on top of the simulated data.
	vec/create x(10) r -0.45 -0.35 -0.25 -0.15 -.05 .05 0.15 0.25 0.35 0.45
* sh  P0
	if [igo]=1 then
	  vec/create y(10) r 78. 106. 161. 305. 468. 468. 305. 161. 106. 78.
* sh  P0'
	elseif [igo]=2 then 
	  vec/create y(10) r 54. 56. 61. 64. 66. 66. 64. 61. 56. 54.
	endif
* sh  Very last number on next line should equal ymax; change if needed
	hi/cre/2dhisto 42 ! 100 -0.5 0.5 100 0. [ylim]
	set hcol 2
	vec/hfill y%x 42
	set mtyp 22
	hi/pl 42  S
	key 0.15 0.9*[ylim] 22 'Real Data'
	set hcol 1
* sh
*
    (shsig)  opt nsta; igset mtyp 1; 
	     title 'Resolution (sigma) vs. Relative Position on Strip, P0 Test Chamber'; 
             nbins = 12
             vec/cre par(3); vec/cre err(3); vec/cre av([nbins]); 
	     vec/cre eav([nbins])
             vec/cre sigmf([nbins]); vec/cre esigmf([nbins])
             id =[ch][dd]
             id2 = [id]-2
             igset mtyp 1
             sliy [id] [nbins]
             hi/proj [id]
             opt stat
             opt fit
             nx  = $HINFO([id],'XBINS')
             xmn = $HINFO([id],'XMIN')
             xmx = $HINFO([id],'XMAX')
             dx  = ([xmx]-[xmn])/[nx];
             sigma xbin = array([nbins],[xmn]+[dx]/2#[xmx]-[dx]/2)
        do ib = 1, [nbins]
          hi/cop [id].sliy.[ib] 99  
*          hi/fit 99 g ! 0 par ? ? ? err
          hi/fit 99 g 0 0 par ? ? ? err
*             rms =   $HINFO(99,'RMS')
*             mean =  $HINFO(99,'MEAN')
*             xxmin = [mean]-2*[rms]
*             xxmax = [mean]+2*[rms]
*          hi/fit 99([xxmin]:[xxmax]) g ! 0 par ? ? ? err
*          wait
          hi/fit 99([xxmin]:[xxmax]) g 0 0 par ? ? ? err
          vect/copy par(2) av([ib])
          vect/copy err(2) eav([ib])
          vect/copy par(3) sigmf([ib])
          vect/copy err(3) esigmf([ib])
        enddo
        igset mtyp 1; opt nsta; opt ndat; 
        igset mtyp 20
        sigma av = av; sigma eav = eav
        sigma sigmf = sigmf; sigma esigmf = esigmf
* sh  Plots average distribution of [id] in [nbins]
*        null  [xmn] [xmx] -100. 100.
*	hpl/err xbin av    ? eav    [nbins] 20 0.25
*        atitle 'relative position on strip' 'Average ([m]m)'
* sh  Plots sigma vs. delx
	if ([ch] = 10)  then 
	  ylim = 250.
	elseif ([ch] = 11)  then 
	  ylim = 250.
	elseif ([ch] > 11)  then 
	  if [igo]=1 then
	    ylim = 500. 
	  else 
	    ylim = 100.
	  endif
	endif
	null  [xmn] [xmx] 0. [ylim]  
        hpl/err xbin sigmf ? esigmf [nbins] 20 0.25
	key -0.45 0.9*[ylim] 20 'Simulation'
        atitle 'relative position on strip' '[s] ([m]m)'
* sh  This plots the experimental data on top of the simulated data.
	vec/create x(10) r -0.45 -0.35 -0.25 -0.15 -.05 .05 0.15 0.25 0.35 0.45
* sh  P0
	if [igo]=1 then
	  vec/create y(10) r 78. 106. 161. 305. 468. 468. 305. 161. 106. 78.
* sh  P0'
	elseif [igo]=2 then 
	  vec/create y(10) r 54. 56. 61. 64. 66. 66. 64. 61. 56. 54.
	endif
* sh  Very last number on next line should equal ymax; change if needed
	hi/cre/2dhisto 42 ! 100 -0.5 0.5 100 0. [ylim]
	set hcol 2
	vec/hfill y%x 42
	set mtyp 22
	hi/pl 42  S
	key 0.15 0.9*[ylim] 22 'Real Data'
	set hcol 1
* sh
*
 (part) zone; title 'Particle type distribution';
                 label 1 6 [g] e^+!   e^-!  [n] [m]^+!   [m]^-! 
                 set ndvx -6.15    
                 set vsiz 0.45
                 set xval 0.6
                 set dmod 1
                 opt logy
                 set *wid 6
                 yy = 1.e6
                 dy = 2
                 null 1 7 1 [yy]
                 c = 1
                 i = 0
                 for id in 312420 300010 300020 300030 300040
                   set hcol [c]
                   hi/pl [id] s
                   yy = [yy]/[dy]
                   if [id]=312420 then
                     txt = 'CSC'
                   else
                     i = [i]+1
                     txt = 'SC-'//[i]
                   endif
                   key 2.0 [yy] [c] [txt] 1.5 F 
                   c = [c]+1
                   if [c]=5 then
                     c = 6
                   endif
                 endfor
                 set ndvx
                 set hcol 1
    (mom) zone; title 'Momentum distributions';
                 zone 1 2
                 set *wid 6
                 set stat 1100
                 set vsiz; set xval
                 null 0 3 1. 1e4
                 set hcol 2; hi/pl 312425 s;  key 2.0 2.e3 2 [m]^+! 1.5 F
                 set hcol 3; hi/pl 312426 s;  key 2.0 1.e3 3 [m]^-! 1.5 F
                 atitle 'P (GeV/c)' ''
                 null 0 0.05 1. 1e4
                 set hcol 2; hi/pl 312423 s;  key 0.02 2.e3 2 e^-! 1.5 F
                 set hcol 3; hi/pl 312422 s;  key 0.02 1.e3 3 e^+! 1.5 F
                 atitle 'P (GeV/c)' ''
endcase
close 67
return