*
*	Kumac file to extract fitting paramaters
*	for the E-counters time spectrae
*	LYM 21-oct-1997
*
opt stat
set stat 111111
opt fit
set fit 111111
zone 2 2
*
* Creation of the various vectors
*	par(4)= 1: background- 2,3,4 gaussian (amplitude,centroid,sigma) 
*	centroids(8,384): output vector
*	min_bound, max_bound: limits of integration defined by sigma.
*	scaler(64): vector in which the raw scalers are dumped
*	scal_bound: number of sigmas about the centroids we'll be
*		integrating over
*	
v/creat par(4)
v/creat centroids(8,384) 
v/creat temp(1024)
v/create test1
v/create test2
v/create test3
v/creat min_bound
v/creat max_bound
v/creat scal_bound
v/creat oth_scal
v/creat integral
v/creat background
v/creat max_val(1024)
v/creat s1
v/creat max_gau
v/creat scaler(64)
v/inp scal_bound 2.5
v/inp oth_scal 4.
*
* Get input data file and go into tag subdirectory
*
read filename 'Please enter rzn filename'
h/file 1 [filename]
title [filename]
cd tag
*
* Read scaler histo first, for normalization purposes
*
h/get/cont 401 scaler
*
* Analysis of all 384 e-counters
*
do i=1,384
*
* centroids(1,x) is the e-counter number
*
v/inp centroids(1,[i]) [i]
*
* starting parameters for fit
*
v/inp par 4. 30. 1035. 10.
h/fit [i](900.:1100.) p0+g 0E 4 par !
*
* testing parameters to see if fit worked
*
v/copy par(3) test1
v/copy par(1) test2
v/copy par(2) test3
*
* Problem finding the peaks. Try again with different
* starting parameters
*
if (test1.le.1000 .OR. test1.ge.1060) then
v/inp par 4. 30. 1040. 10.
h/fit [i](900.:1100.) p0+g 0E 4 par !
endif
if (test2.lt.0 .OR. test3.lt.0) then
v/inp par 4. 30. 1040. 10.
h/fit [i](900.:1100.) p0+g 0E 4 par !
endif
*
* Integrate over peak +/- (scal_bound*sigma)
*
if (par(4).gt.0) then
v/op/vmul par(4) scal_bound min_bound
v/op/vadd par(3) min_bound max_bound
v/op/vsub par(3) min_bound min_bound
else
v/op/vmul par(4) scal_bound min_bound
v/op/vsub par(3) min_bound max_bound
v/op/vadd par(3) min_bound min_bound
endif
*
* dump values of histogram into vector temp
*
h/get/con [i] temp
*
* Initialized integral and background variables
*
v/inp integral 0
v/inp background 0
*
* determine vector location as a function of tdc value
*
v/op/vdiv min_bound oth_scal i_min
v/op/vdiv max_bound oth_scal i_max
if (i_min.lt.1) then
v/inp i_min 1
v/inp i_max 1023
endif
if (i_max.gt.1023) then
v/inp i_min 1
v/inp i_max 1023
endif
sigma i_min=int(i_min)
sigma i_max=int(i_max)
*
* perform the integral
*
do k=i_min,i_max+1
v/op/vadd integral temp([k]) integral
v/op/vadd background par(1) background 
enddo
*
* subtract the integral from the noise = get real signal
*
*v/op/vsub integral background integral
*
* Load various value into output vectors
*
h/get/conten [i] max_val
sigma s1=vmax(max_val)
v/op/vsub s1 par(1) max_gau
v/inp centroids(2,[i]) par(1)
v/inp centroids(3,[i]) max_gau
v/inp centroids(4,[i]) par(3)
v/inp centroids(5,[i]) par(4)
v/inp centroids(6,[i]) integral
v/inp centroids(7,[i]) background
v/inp centroids(8,[i]) scaler(64)
h/plot [i](900.:1100.)
enddo
*
* dump output to file
*
v/write centroids centroids.dat 8(1x,g12.5)
*v/write centroids centroids.dat '7(1x,f8.2),1x,g12.4'