;+
; NAME:
;	BCS_PHA
; PURPOSE:
;	Extract and plot PHA data from BCS Radioactive Source Calibration
; CALLING SEQUENCE:  
;
; INPUT:
;
; OUTPUT:
;
; OPTIONAL OUTPUT:
;
; HISTORY
;	RDB, 1989, with enhancements
;	16-Feb-92	Asks for approx. start time
;-
chan_id  = ['Fe XXVI','Fe XXV ','Ca XIX ','S XV   ']
spec_id  = ['BCS-A','BCS-B']
temp_ahk = [95,96]		;pointers for BCS crystal temps conversion
hv_ahk   = [32,33]		;pointers for BCS crystal HV's conversion

ps_reset
!p.multi=0

ansdf='N'
;read,'* Use defaults [def=Y] ',ansdf
;ansdf=strupcase(strmid(ansdf,0,1))
;if ansdf ne 'N' then ansdf='Y'

ans = ''

break_file,infil(0),disk_log,dir,filnam,ext,fver	;????

;;tarr = intarr(7)
;;print,' '
;;read,'* Enter approximate start time [h,m,s,ms,d,m,y]: ',tarr
;;print,'This corresponds to record: ',tim2dset(dp_sync,tarr)

;	cal. source cal. when Cal. B enabled
caben = (dp_sync.specstatus and 1)
;;;qq = where(caben ne 0)
;;;ista = (min(qq)/4+1)*4
;;;istp  = max(qq)/4*4
cabenx = where(caben(1:*)-caben(0:*) ne 0)
print,' '
print,'CALB enabled between -'
print,cabenx/2
ista=0 & istp=0
read,'* Enter region to be examined: ',ista,istp
print,' '

;;;;ista = ista*2 & istp = istp*2
fi = dp_sync.dp_time(3)
ista = min(where(fi(ista:*) eq 0)) + ista
nnpha = (istp - ista)/4
pha = byte(dp_sync.pha_counts,ista*8,32,nnpha)		;8 PHA bytes/rec
help,pha

;	check if RBM on  (RBM is once/frame)
rbm = (dp_sync.rbm_status and 128)/128
rbmx = where(rbm gt 0)/4

;	check for the channel changes
phachans = dp_sync.pha_control/64
phcnx = all_vals(where(((phachans(1:*) - phachans(0:*)) ne 0))/2)
phcnx = [ista,phcnx(where((phcnx ge ista) and (phcnx le istp)))]

if phcnx(0) eq -1 then phcnx = [0,n_elements(phachans)-2]
;;if phcnx(0) gt 0 then phcnx = [0,phcnx]
print,phcnx

cstart=0 & cstop=3
;;istart = min(where(bit_map ne 255)) 
;;istop = max(where(bit_map ne 255))
istart = 0 & istop = n_elements(pha)/4+ista	;??

if ansdf eq 'Y' then begin
cchans = phachans(istart)	;!!!!
cstart = cchans & cstop = cchans
	endif	;defaults

for jchn = cstart,cstop do begin
  vt200
;;  chx = where(phachans(8:istop/8*8-1) eq jchn)
  chx = all_vals(where(phachans(ista*2:istp*2) eq jchn)/2) + ista
  stapos = min(chx) + 3*4 	;;sync?????    ((min(chx)/8)+2)*8
  stppos = min([istop,phcnx(min(where(phcnx ge stapos+1) ))])
  print,stapos,stppos
  starec = stapos
	istart = starec
  stprec = stppos -1   
	inum = (stprec - starec)/4 -1 ;???
  print,starec,stprec,istart,inum

if inum le 0 then goto,nextch

try_again:
;;print,'Istart,inum',istart,inum
print,' '
print,'The numbers below show when the PHA channel changed'
print,'Because of telemetry dropouts, etc, there may be false changes'
print,phcnx
print,' '
print,'Data for channel',jchn,'  starts at frame:',istart
print,'Number of PHA frames before channel change is:',inum
print,'To allow for telemetry delays, it is best to choose new limits'

ans=''

if ansdf ne 'Y' then begin
ii1=0 & ii2=0
read,'Enter modified values (-1,-1 for above) ',ii1,ii2
if ii1 gt 0 then begin
	istart=ii1 & inum=ii2
	starec = istart   & stprec = starec + inum
	stapos = starec & stppos = (stprec+1)
	endif

;	flag RBM activity!
if n_elements(rbmx) gt 1 then begin
;;	print,starec,stprec,min(rbmx),max(rbmx)
if ((min(rbmx) ge starec) and (min(rbmx) le stprec)) $
	or ((max(rbmx) ge starec) and (max(rbmx) le stprec)) then begin
	print,' '
	print,'**************'
	print,'*** RBM active in this interval between',min(rbmx),max(rbmx)
	print,'**************'
	goto, try_again
	endif
	endif

;	get main title????
main_title = 'BCS PHA Data'
mtitle = main_title
;;mtitle=''
;;lent = strlen(main_title)
;;if lent gt 0 then begin
;;	print,'Main title is: ',main_title
;;	read,'Do you wish to change this ',ans
;;	ans = strupcase(strmid(ans,0,1))
;;	mtitle = main_title
;;	if ans eq 'Y' then goto, getmain
;;	endif
;;goto, gotmain	
;;getmain:
;;read,'* Enter Title ',mtitle
;;gotmain:
	endif	;defaults
print,' '

;	get start time
;;idx2list,dp_sync(istart),['time','day'],fxtime
int2ex,dp_sync(istart).time,dp_sync(istart).day,tarr
fxtime = fmt_tim(tarr)
print,fxtime

;	get PHA channel
k = all_vals(dp_sync(stapos:stppos-1).pha_control/64)
if n_elements(k) ne 1 then begin
	print,'Error PHACH',k
	read,'Continue',ans
	endif
phach = k(0)
ispec = phach/2		;spectrometer A or B

;	get PHA time
k = all_vals(dp_sync(stapos:stppos-1).pha_control and 63)
if n_elements(k) ne 1 then begin
	print,'Error PHATM',k
	read,'Continue',ans
	endif
phatm = k(0)

print,'PHACH, PHATM',phach,phatm
print,spec_id(ispec),'   ',chan_id(phach)

;	get HVPS trim value (depending on PHA channel)
jidx = sel_qs(qs.instr,dp_sync(stapos:stppos))
k = all_vals(jidx)
if n_elements(k) ne 1 then begin
	print,'Error HVTRIM, SEL_QS index changes',k
	ansz=''
	read,'Continue?:',ansz
	endif
trims = qs.instr(jidx(0)).hv_control 
trims = trims(0)
hvtrim = trims and 7
if ispec eq 1 then hvtrim = trims/16 and 7

;	get HVPS voltage value of this unit
ehtval = ahk_conv(dp_sync.hv_mon(ispec),hv_ahk(ispec),fi)	;,dp_sync)
ehtval = ehtval(stapos/2:stppos/2-1)
k = all_vals(ehtval)
if k(0) eq 0 then k = k(1:*)
if n_elements(k) ne 1 then begin
  print,'Error EEHTX',k
  endif
ehtval = fix(k(0))

;	get temperature of this unit
tempr = ahk_conv(dp_sync.crytemp(ispec),temp_ahk(ispec),fi)	;,dp_sync)
tempr = tempr(stapos/2:stppos/2-1)
tempr = total(tempr)/n_elements(tempr)

Print,'HV Trim and mon, Cryst temp.',hvtrim,ehtval,tempr

istartxx=(istart-ista)/4
phas = pha(*,istartxx)
;;if inum gt 1 then phas = sum(pha(*,istartxx:istartxx+inum-1),1)/float(inum)
if inum gt 1 then begin
	zz = pha(*,istartxx:istartxx+inum-1)
	cc = zz(*,0)*0		;ss = size(zz)  is ss(2) == inum??
	for c=0,inum-1 do cc = cc+zz(*,c)
	phas = cc/float(inum)
	endif
xx=indgen(32)

;	do not fit end bins
;;plo = 5 & phi = 29
plo = 8 & phi = 30
;;sd = xstdev(phas(5:29),xb) & xb=xb+5
;;xw = 3. & plo = fix(xb)-2*xw & phi = plo+4*xw	;xw=sigma, use 2 sig.
yy = gaussfit(xx(plo:phi),phas(plo:phi),a);,0)	;get fit parameters
print,a
;;;	narrow the fit to 7* half width
;;plo = a(1) - 7*a(2) & phi = a(1) + 7*a(2)
;;yy = gaussfit(xx(plo:phi),phas(plo:phi),a);,0)	;get fit parameters
;;print,a

fwhm = 2.*sqrt(2.*alog(2.))*a(2)
print,'fwhm',fwhm,fwhm/a(1)*100.
;;xw = fwhm & plo = a(1)-2*xw & phi = a(1)+2*xw
;;yy = gaussfit(xx(plo:phi),phas(plo:phi),a);,0)	;get fit parameters
;;print,a
;;fwhm = 2.*sqrt(2.*alog(2.))*a(2)

if !d.name eq 'TEK' then  read,ans		;???????????

TEK
!p.multi = 0		;single plot
;;reset_plot,!device
!fancy=1
!mtitle = mtitle
;;!mtitle = 'BCS PHA Calibration using: Fe 55'	;K!4a!3'
!ytitle = 'sum(pha(*,ista:ista+inum-1),1)/inum
!xtitle = 'PHA bin'
set_xy

;;xx1=findgen(250)/10.+5
xx1=findgen((phi-plo)*10)/10. + plo
funct,xx1,a,yy1 	;generate curve
!xticks=8 ;;;& !type=4
set_xy,0,32,0,max([max(yy1),max(phas(0:29))]) & !c=0
!psym=0 & plot,xx1,yy1,xstyle=1
;;;!psym=10 & oplot,xx,phas,xstyle=1
!psym=10 & oplot,xx,phas
!type=0 & !xticks=0

;	create some info
ident = '  ('+ spec_id(phach/2) + ', ' + Strtrim(chan_id(phach))+')'

;?????? needs fixing
;;durn = (s_tim(istart+inum) - s_tim(istart))*3600.
durn = (phatm*0.5-0.125)

set_xy,0,1,0,1
;
xyouts,.7,.95 $
	,string('FWHM =',fwhm,fwhm/a(1)*100.,'(a,f7.3,2x,1h(,f5.2,2h%))') $
	,siz=1.
xyouts,.7,.91,string('Centroid =',a(1),'(a,f7.3)'),siz=1.
;
xyouts,.03,.95,string('PHA Channel =',PHACH,ident,'(a,i2,a)'),siz=1.
xyouts,.03,.91 $
	,string('HV Trim =',HVTRIM,'  HV =',EHTVAL,' V','(a,i2,a,i5,a)') $
	,siz=1.
;;frxtime = frtime(*,istart*8)
;;frxtime = intarr(6)			;************
;;mm = dp_sync(istart).time
;;frxtime(3) = mm(0:2)
;;;xyouts,.05,.87,string('Start Time =',frxtime(3:5),'(a,i3,2(1h:,i2))')
;;xyouts,.03,.87,string('Start Time =',frxtime(3:5),frxtime(0:2) $
;;        ,'(a,i3.2,2(1h:,i2.2),1h,,i3.2,2(1h/,i2.2))') $
;;	,siz=1.
;idx2list,dp_sync(istart),['time','day'],fxtime
xyouts,.03,.87,'Start Time = '+fxtime,siz=1.
xyouts,.03,.83,string('Duration =',durn,' secs.','(a,f6.1,a)'),siz=1.
xyouts,.03,.79,string('ISTA,INUM= ',istart,inum,'(a,2i4)'),siz=1.
xyouts,.03,.75,string('Temp =',tempr,' C','(a,f7.2,a)'),siz=1.

;;set_plot,!device
;;set_xy,0,1,0,1
!fancy=-1
xyouts,.76,0.97,filnam+ext,siz=1.2,/normal

nextch:
if ansdf eq 'Y' then goto,endit
if !d.name eq 'TEK' then  read,'Next: ',ans
  if strupcase(strmid(ans,0,1)) eq 'A' then goto,try_again
  if strupcase(strmid(ans,0,1)) eq 'E' then goto,endit

  endfor	;loop through channels

endit:
set_xy
end