131 lines
2.5 KiB
Plaintext
131 lines
2.5 KiB
Plaintext
*
|
|
* $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
|
|
*
|
|
*-------------------------------------------------------------------------------
|
|
*
|
|
|
|
|
|
|