pro readcal,filename,xm,ym,nne,eav,gausspeak
;
;	corrections:	long integer ii 6/3/96
;	attempt to get outfile in right order 6/3/96
;
;	attempt to 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
;	read all this data
; 	
;	Print header info 1/29/97
;
;		Version 1.1	adding looping over long data files,
;				build up histograms of positions, energy
;				get average energy, number of pedestals,
;				etc.

;	BUG: last time through loop it overshoots the file.  Not a showstopper.
;	use for floodtests and high voltage files


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

;	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)		; data arranged in groups of 4
ii=0L
;readu,3,data
repeat begin
	readu,3,count
	ii=ii+4
endrep until eof(3)
			; counted how many events in data file

;data=bytarr(ii)
close,3
openr,3,filename	; read header again and ignore
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
;
;		Since files are too large to easily manipulate, loop
;		through subset of the file, and calculate histograms
;		of parameters of interest	1/30/97

l_unit=100000	;analyze l_unit events at once

data=bytarr(l_unit)	;pick a chunk

l_end=ii/l_unit +1	; find out where to end loop

histall=fltarr(641,641)

for i=0, l_end do begin		;META-loop begins *******************

start=l_unit*i
last=l_unit*(i+1)-1
if(last gt ii) then last=ii-1
print,'current chunk covers: ',start,last

readu,3,data	;read a chunk of the data file



;n=ii/8		;old way, when all data was being looked at simultaneously

n=(last-start+1)/8


;	For flood test data, the arrays quickly get too big
;	to play fast and loose with copying them.  Break the
;	data array up into subsets at this point, and gather
;	statistics on each chunk.
;					PTB	1/27/97



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

data1=reform(data,2,4*n)

data1=long(data1)		;RUNS OUT OF MEMORY HERE!

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
;	apparently 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

if total(locp) ne -1 then peds=xy(*,locp) else peds=0

if total(locxy) ne -1 then ev=xy(*,locxy)

;	From ADC values of pedestal events, get four offsets

ptot=lonarr(4)
if peds eq 0 then ptot(0:3) =0
if peds ne 0 then for j=0,3 do ptot(j)=total(peds(j,*))
if total(ptot) ne 0 then ptot=ptot/n_elements(locp)

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

;	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,*))



hx=histogram(xm,min=0,max=640)

hy=histogram(ym,min=0,max=640)

histo,xm,ym,histout


histall=histall+histout

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

 e=(ev(0,*)+ev(1,*)+ev(2,*)+ev(3,*))/32
;e=total(ev)/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
tvscl,histall
endfor			;META loop ends *********************

close,3

end
@histo