*
*	$mcv3k_root\paw\plot_range.kumac
*
*	Plot mu+ ranges in Al up to p_in = 29.8 MeV/c and fit
*	a1xp**3.5 to the data.
*	Ranges calculated with mcv3k. 
*	Histograms saved in $mcv3_root\data\al_range.rz.
*
*	requires one input parameter:
*
*	[1]	material, for example al
*
*	TP,	28-Nov-2000,	PSI LEM group	
*
*-----------------------------------------------------------------------------
*
if ( [#] = 0 ) then
	mess
	mess Need input parameter to build input file name.
	mess For example:
	mess
	mess plot_range al
	mess 
	mess Input file name will be $mcv3k_root/data/al_range.rz
	mess
	exitm
endif
*
mat = [1]
*
*------------------------------------------------------------------------------
*
close 0 | close all open units
hi/de 0 | delete all histograms in //pawc
v/de *  | delete all vectors
*
v/crea hid(100)	| vector for histogram ID's
v/crea x0(2)    | for defining graphic window
v/crea y0(2)	| same for y
*
*-------------------------------
*
file = $mcv3k_root/data/[mat]_range.rz | file with range histograms
if ( $fexist([file]) = 0) then
	mess 
	mess File [file] does not exist ?!
	mess
	exitm
endif
*
*	open file
*
hi/file 20 [file]
hrin 0
close 20
nhist = 0
do i = 1, 300	| look for histograms; ID = 10 x momentum [MeV/c]
	if ( $hexist([i]) ) then
		mess Found histogram [i]
		nhist = [nhist] + 1
		v/inp hid([nhist]) [i]
	endif
enddo
*
* now we are ready to extract information from histograms
*
v/crea r([nhist])	| vector for mean ranges
v/crea er([nhist])	| vector for RMS range
v/crea p([nhist])	| vector for momenta
v/crea ep([nhist]) r 8*0. | error in momentum = 0
*
do i = 1, [nhist]
	id = hid([i])
	mom  = [id]/10.
	range = $hinfo([id],'MEAN')
	rms   = $hinfo([id],'RMS')
	v/inp r([i]) [range]
	v/inp er([i]) [rms]
	v/inp p([i]) [mom]
enddo
xmax = [mom]*1.15
ymax = [range]*1.07
v/inp x0(1) 0.
v/inp y0(1) 0.
v/inp x0(2) [xmax]
v/inp y0(2) [ymax]
*
*	plot
*
set ysiz
zone
opt grid
mat = $upper([mat])
title  [mat]//', [m]^+! mean range, [D]p=0'
*
*
*
graph 2 x0 y0 aw
hpl/err p r ep er [nhist] 20 0.3
*
*	fit to data par(1)xp**3.5
*
v/crea par(1) r 1.
opt nstat
opt fit
set fit 111
v/fit p r er $mcv3k_root/paw/range.f s 1 par
scale = par(1)
max = p(1)
func/pl [scale]*x**3.5 0. [max] s
atitle 'p (MeV/c)' 'mean range ([m]m)'
*
xpos = 0.05*x0(2)
ypos = 0.95*r([nhist])
str = $mcv3k_root/paw/plot"-#range.kumac
text [xpos] [ypos] [str] 0.2
ypos = 0.90*r([nhist])
text [xpos] [ypos] 'Fit function "J# P1 x p^3.5!' .25
*
*
*
title ' '
opt ngri
opt stat
*
*-------------------------------------------------------------------------------
*