INTRODUCTION

# startup
sm  - to enter sm
quit - to exit sm
help - get help in sm
apropo - find related topics

# define a vector like this, access as rr o rr[1] rr[2] etc
# also can use set vec=start,end,delta to generate a series
set rr = random(10)       
print{rr}                 # will print values
# define a scalar like this
define tt   3.14 
# read in data from file.dat like this
data "./file.dat"
read {u 1 v 2 w 3}  # u from col 1 etc.

# set up a basic plot
device x11 # o device postencap "filename.ps"
erase
limits xmin xmax ymin ymax # o limits x y    if x,y are vectors to plot
box # draw box around plot

# actual plotting
connect xvec yvec # join points
points xvec yvec # plot points
# histograms
histogram xvec yvec # Makes a histogram of height yvec at xvec.
# quick errors
error_x, error_y xvec yvec evec
# error bars
errorbar o logerr


# random numbers: poisson errors
define size 20 # size of sample
# later repeat for size=200 points - whats better?
set rr=random($size)
set bins=0,1,0.1   # 10 bins between 0 and 1
set rrhist=histogram(rr:bins)
erase
lim bins 0 10
limits bins rrhist 
box
xlabel Value of random 
ylabel No. samplesç
# plot muestra
histogram bins rrhist
angle 45 shade histogram 500 bins rrhist
angle 0
# now add error bars : poisson
set rrhister=rrhist**0.5  # root error bars
error_y bins rrhist rrhister

# now add madre distribucion: for size numbers, we expect size/10 per bin
set eey=rrhist*0 + $size/10.
ltype 1 connect bins eey  # expectation in each bin
ltype 0
# now add error bars for madre : poisson
set eeyer=eey**0.5  # root error bars
#error_y (bins+0.01) eey eeyer # lets offset to make it clear

# now calculate Chi-square!
# start either from the parent distribution or the muestra
set mean=  {0.5} # this is the mean of madre
set mean = sum(rr) / $size # this is the mean of the muestra
print{mean}
set diff = (rrhist - mean)**2  # difference square
set diff2 = diff / (rrhist+0.0001)  # in case y=0 in any bin
set chisq = sum(diff)
# reduced chisq: divide by deg of freedom (n. of points (-1 if you used mean))
set chisqred=chisq/$size
print{chisqred}

# compare this to the tables of chisq


# lets instead compare to a gaussian
# take mean and sigma from the muestral data
set gbins= 0,1,0.01
set mean = sum(rr) / $size # this is the mean of the muestra
set diff2 = (rr - mean)**2 # difference square
set sig = (sum(diff2) / $size)**0.5
# the peak is normalized to the width and mean (total area = 1)
set gaussy= (1.0/(sig*(2*3.14)**0.5)) * exp( (-0.5 * ((gbins - mean)/sig)**2))
# overplot gaussy for comparison
connect gbins (gaussy*$size*0.1)

# reduced chisq for gaussian comparison
set gaussy= (1.0/(sig*(2*3.14)**0.5)) * exp( (-0.5 * ((bins - mean)/sig)**2))
set gaussy= gaussy*$size*0.1
set diff = (rrhist - gaussy)**2  # difference square
set diff2 = diff / (rrhist+0.0001)
set chisq = sum(diff2)
# reduced chisq: divide by deg of freedom (n. of points (-1 if you used mean))
set chisqred=chisq/$size
print{chisqred}

# compare this to the tables of chisq