pro readcal1,p,q,filename,xm,ym,nne,eav,gausspeak
;
;	HAZARD: This was written VERY Quickly and may change drastically!
;	VERSION: 0.0					PTB   2/3/97
;	----------------------------------------------------------------
;	corrections:	long integer ii 6/3/96
;	get outfile in right order 6/3/96
;
;	define a region of interest around each point
;	p and q are coordiates on the detector

;	Choose 40 as ROI to fit Kevin's program

;	use to analyse an entire directory of positionertest data
;
;	NOTE: use this and FILEREAD together, ONLY on positioner files

;	-----------------------------------------------------------------


;Gather and print header info

head=string(333)
openr,3,filename
readf,3,head
print,head
readf,3,head
print,head
readf,3,head
print,head
readf,3,head
print,head
readf,3,head
print,head
readf,3,head
print,head
readf,3,head
print,head
readf,3,head
print,head
readf,3,head
print,head
readf,3,head
print,head
readf,3,head

print,'end of header info
stop

;	Now, the first chunk of events are pedastal events, and there is
;	a LAST_PED marker (1000).  The 16 bit events are formatted as in
;	Kevin's EVENT PROCESSING document.  BITS 0-11 are the ADC values,
;	and BITS 12-15 contain various flags in the four words (i.e.,
;	anti-co, risetime events, upperlevel, and pedastal events.
;	*****************************************************************
data=bytarr(20000)
count=bytarr(4);	count up units of 4 ADC values in file (cumbersome)
ii=0L
;readu,3,data
repeat begin
	readu,3,count
	ii=ii+4
endrep until eof(3)


data=bytarr(ii)
close,3
openr,3,filename
readf,3,head
readf,3,head
readf,3,head
readf,3,head
readf,3,head
readf,3,head
readf,3,head
readf,3,head
readf,3,head
readf,3,head
readf,3,head

readu,3,data
close,3



locdat=where(data gt 0)

close,3
n=ii/8

data1=data
 byteorder,data1  ; for now, assume all must be swapped.

data1=reform(data1,2,4*n)
data1=long(data1)
new=data1(0,*)*256+data1(1,*)	; these are now the 16-bit words 
new1=reform(new,4,n)		; now in order as expected

;	calculate x and y position from xy array
;	organized as [y-,y+,x-,x+]

xy=new1 and 4095	;bitwise AND takes off last 12 bits
head=new1/4096		;this leaves head info only


;	If the first of the four header bytes is 1101 (13) then this
;	is a pedestal event. Else it is 1100 (12), which is a science
;	event.
;	***************************************************************

locp=where(head(0,*) eq 13)
locxy=where(head(0,*) eq 12)
locbad=where(head(0,*) ne 12 and head(0,*) ne 13)


;	Now pedestal events are averaged to get the DC offsets to subtract
;	from normal events

peds=xy(*,locp)
ev=xy(*,locxy)

;	From ADC values of pedestal events, get four offsets

ptot=lonarr(4)
for i=0,3 do ptot(i)=total(peds(i,*))
ptot=ptot/n_elements(locp)

for i=0,3 do ev(i,*)=ev(i,*)-ptot(i)

;	Using x(+) or y(+) instead of the difference in the denominator
;	simply replaces a scale from (-1,1) with one from (0,1)
;

ym=ev(1,*)*640.
ym=ym/(ev(1,*)+ev(0,*))

xm=ev(3,*)*640.
xm=xm/(ev(3,*)+ev(2,*))

print,'p and q are',p,q

lxy=where(xm lt p+40 and xm gt p-40 and ym lt q+40 and ym gt q-40)
;ly=where(ym lt q+40 and ym gt q-40)

xuse=xm(lxy)
yuse=xm(lxy)

hx=histogram(xuse,min=0,max=640)
hy=histogram(yuse,min=0,max=640)



;	calculate energy as the sum of the x and y over 32

e=(ev(0,lxy)+ev(1,lxy)+ev(2,lxy)+ev(3,lxy))/32

loce=where(e ge 0 and e le 100000)	;choose energy region
nne=n_elements(loce)
esum=total(e(loce))
eav=float(esum)/float(nne)
egauss=e(loce)

;	Get set to fit the energy distribution with Gaussian

yax=histogram(egauss)
xax=findgen(n_elements(yax))

result=gaussfit(xax,yax,param)
gaussp=max(result,gausspeak)

;print,nne,eav,gausspeak

end

;	known bugs: getting an underflow to the gauss fit once in a while