;+
; NAME:
;       ANAL_BCS_PHA
; PURPOSE:
;       Extracts BCS PHA data and analyses it
; CATEGORY:
; CALLING SEQUENCE:
;	anal_bcs_pha,fileid,dp_sync,qs
; INPUTS:
;       fileid          ID of file being processed
;       dp_sync         dp_sync data of this file
;       qs              qs data from this file
; OPTIONAL INPUT PARAMETERS:
; KEYWORD PARAMETERS:
;       plot            if set, plots result to screen
;                        if !d.name='PS', plots to "phFILNAM.ps" on curr. dir.
;       log             if set, writes log to "gain_log.dat" on curr. dir.
;       verbose         if set, print aditional information
; OUTPUTS:
; OPTIONAL OUTPUT PARAMETERS:
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
;	Only tested with source cals
; PROCEDURE:
;	Finds Source Cal within suppilied data, extracts an analyses
;	Logs results if required
; MODIFICATION HISTORY:
;	RDB  26-Apr-94	Written
;	RDB  19-May-94	Documented; return (after logging) in no cal.
;	RDB  22-May-94  Fixed bug that stopped remote execution
;	RDB   1-Jun-94	Trapped cases (tests?) where only single channel
;       RDB   3-Aug-94  Only use PHA bins 0-30 when determining max. value
;	RDB   9-Feb-95	Pick up log name from CAL_LOG if defined...
;-

pro anal_bcs_pha,fileid,dp_sync,qs,plot=plot,verbose=verbose,log=log

if keyword_set(plot) then begin
  if !d.name eq 'X' then loadct,4
  clearplot
endif

print,fileid
help,dp_sync,qs.instr

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

fi = dp_sync.dp_time(3)

cpu0 = where(dp_sync.dp_time(3) eq 0)
ncpu = n_elements(cpu0)         ;/4*4           ;??? whole sets of dp_sync
print,cpu0(0)

cpu0_time = anytim2ints(dp_sync(cpu0))

pha = bytarr(32,ncpu)
for j=0,ncpu-2 do $
  for jj = 0,3 do pha(jj*8,j) = dp_sync(cpu0(j)+jj).pha_counts


pha_cnts = total(pha,1)/.375
pha_chn = dp_sync(cpu0).pha_control/64

help,pha,cpu0_time,pha_cnts,pha_chn

caben = (dp_sync(cpu0).specstatus and 1)
cabenx = where(caben(1:*) -caben(0:*) ne 0)/2
;	may be a problem if cal not found.
if n_elements(cabenx) ne 2 then begin
  print,'** This file does not contain Source Calibration data **  ',fileid
  if keyword_set(log) then begin
    openw,lune,'gain_log.err',/get_lun,/append
    printf,lune,'** Problem with file, no source cal data:  ',fileid
    free_lun,lune
  endif
  return
endif

print,' '
print,'CALB enabled between -',cabenx,format='(a,2i6)'
w_cal = indgen(cabenx(1)-cabenx(0)+1)+cabenx(0)

sdate = anytim2ex(cpu0_time(w_cal(0)))

;	determine times of channel changes, and remove funnies...
channels = pha_chn(0,w_cal)
print,'Channel:',all_vals(channels)
ccc = channels(1:*)-channels(0:*)
chan_c = [cabenx(0),where(ccc ne 0)+1+w_cal(0),cabenx(1)]

for kkk = 0,3 do begin
  if chan_c(1)-chan_c(0) lt 4 then chan_c = chan_c(1:*)
  ncc = n_elements(chan_c)
  if chan_c(ncc-1)-chan_c(ncc-2) lt 4 then chan_c = chan_c(0:ncc-2)
endfor

print,chan_c

;	search for the channel limits
if keyword_set(verbose) then plot,w_cal,(total(pha,1))(w_cal)

limits = intarr(2,4)-1		;-1 value in case channel not found

for jj=0,n_elements(chan_c)-2 do begin
  lims = [chan_c(jj)+3,chan_c(jj+1)-1]
  ccx = all_vals(pha_chn(0,lims(0):lims(1)))
  if keyword_set(verbose) then begin
    oplot,[lims(0),lims(0)],[0,800]
    oplot,[lims(1),lims(1)],[0,800]
    print,lims,ccx
  endif
  ccx = ccx(0)
  clim = limits(1,ccx)-limits(0,ccx)+1
  nlim = lims(1)-lims(0)+1
;  print,clim,nlim
  if nlim gt 8 and nlim gt clim then limits(*,ccx)=lims
endfor

print,limits

;	analyse the channels with the selected limits
ps_file = 'ph' + strmid(fileid,3,6) + '.ps'
print,'PS_FILE= ',ps_file
if keyword_set(plot) and !d.name eq 'PS' then begin
  device,filename=ps_file
  device,/land
endif

!p.multi=[0,2,2]

results = fltarr(2,4)-1.	;-1 for skip when plot
temps = fltarr(4)
xx = indgen(32)

for jchan = 0,3 do begin
 ispec = jchan/2
 ident = '  ('+ spec_id(ispec) + ', ' + Strtrim(chan_id(jchan))+')'
 mtit=string('BCS Channel',jchan+1,format='(a,i2)')

 lims=limits(*,jchan)
 n_spec = lims(1)-lims(0)+1
 print,jchan,lims,n_spec,format='(4i8)'
 cymax = 0
 if n_spec gt 1 then begin
   pha_spec = total(pha(*,lims(0):lims(1)),2)/n_spec
   cymax = max(pha_spec(0:30))		;don't use last bin in max
 endif

 if n_spec eq 1 or cymax eq 0 then begin
  print,'** No Data **',jchan,n_spec,cymax
  if keyword_set(plot) then begin
    plot,xx,xx*0,/nodata,tit=mtit,xra=[0,31],xsty=1,xtit='PHA Bins'
    xyouts,10,.5,'NO DATA',siz=1
  endif

 endif else begin
  if keyword_set(plot) then plot,xx,pha_spec,psym=10, $
    xra=[0,31],xsty=1,xtit='PHA Bins',yra=[0,1.25*cymax],ysty=1, $
    tit=mtit
  w_peak = (where(max(pha_spec(0:30)) eq pha_spec))(0)
  plo = w_peak-6 & phi = w_peak+6
  plo = 8 & phi=30			;!!!!!!!!!!!!!!!!!!!!
  print,w_peak,plo,phi

;	fit the spectrum
  yy = gaussfit(xx(plo:phi),pha_spec(plo:phi),a)
  fwhm = 2.*sqrt(2.*alog(2.))*a(2)
  resolution = fwhm/a(1)*100.
  print,a(1),fwhm,resolution,format='(3f9.2)'
  results(0,jchan) = a(1)		;save results
  results(1,jchan) = fwhm
; 	generate and overplot fitted spectrum
  xx1 = findgen((phi-plo)*10)/10.+plo
  funct,xx1,a,yy1
  if keyword_set(plot) then oplot,xx1,yy1

;	get some needed telemetry points

;	get start time and duration of this channel
  statim = fmt_tim(cpu0_time(lims(0)))
  durn = fix(addtime(cpu0_time(lims(1)+1),diff=cpu0_time(lims(0)))*60+0.5)
  print,statim,durn

;       get temperature of this spectrometer
  tempr = ahk_conv(dp_sync.crytemp(ispec),temp_ahk(ispec),fi)
  tempr = tempr(lims(0)*2:lims(1)*2)
  tempr = total(tempr)/n_elements(tempr)
  print,'Temperature:',tempr
  temps(jchan) = tempr

;       get HVPS trim value of this spectrometer
  jidx = sel_qs(qs.instr,cpu0_time(lims(0):lims(1)))
  k = all_vals(jidx)
  trims = qs.instr(jidx(0)).hv_control
  trims = trims(0)
  hvtrim = trims and 7
  if ispec eq 1 then hvtrim = trims/16 and 7
  print,'Trims:',hvtrim

;       get HVPS voltage value of this spectrometer
  ehtval = ahk_conv(dp_sync.hv_mon(ispec),hv_ahk(ispec),fi)
  ehtval = ehtval(lims(0)*2:lims(1)*2)
  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 = total(ehtval)/n_elements(ehtval)
  print,'High Voltage:',ehtval

;	put in annotation
  if keyword_set(plot) then begin
  set_xy,0,1,0,1
;
  csiz=0.5
  xyouts,.60,.94,string('Centroid =',a(1),'(a,f7.3)'),siz=csiz*1.3
  xyouts,.60,.89 $
        ,string('FWHM =',fwhm,resolution,'(a,f7.3,2x,1h(,f5.2,2h%))') $
        ,siz=csiz*1.3
;
  xyouts,.03,.95,string('PHA Channel =',jchan,ident,'(a,i2,a)'),siz=csiz
  xyouts,.03,.91 $
        ,string('HV Trim =',HVTRIM,',  HV =',EHTVAL,' V','(a,i2,a,i5,a)') $
        ,siz=csiz

  xyouts,.03,.87,'Start Time = '+statim,siz=csiz
  xyouts,.03,.83,string('Duration =',durn,' secs.','(a,f6.1,a)'),siz=csiz
  xyouts,.03,.79,string('ISTA,INUM= ',lims(0),n_spec,'(a,2i4)'),siz=csiz
  xyouts,.03,.75,string('Temp =',tempr,' C','(a,f7.2,a)'),siz=csiz
  if jchan eq 2 then  xyouts,-0.1,-0.25,fileid
  endif
 endelse
  
endfor

!p.multi=0
if keyword_set(plot) and !d.name eq 'PS' then device,/close

if keyword_set(log) then begin
  tempx = temps(where(temps gt 0))
  tempx = total(tempx)/n_elements(tempx)	;average temperature
  report = string(sdate(4:6),'V',results,tempx,fileid,sdate(0:1), $
	format='(i2.2,2(1h/,i2.2),2x,a,x,2f7.3,3(2x,2f7.3),f6.1,2x,a,x,1h(,i2.2,1h:,i2.2,1h))')
  print,report
  if n_elements(where(results lt 0)) gt 4 then begin
    openw,lune,'gain_log.err',/get_lun,/append
    print,'** Problem with file, too few channels: ',fileid
    printf,lune,'** Problem with file, too few channels: ',fileid
    printf,lune,report
    free_lun,lune
  endif else begin
    log_name = get_logenv('CAL_LOG')
    if log_name eq '' then log_name = 'gain_log.dat'
    openw,lunw,log_name,/get_lun,/append
    printf,lunw,report
    free_lun,lunw
  endelse
endif

end