; $Id: dcs_cryst_align.pro,v 1.1 2008/05/20 14:48:01 ymqiu Exp $
;#######################################################################
;
; NAME:
;  dcs_cryst_align
;
; PURPOSE:
;  This program plots crystal alignment files
;
; CATEGORY:
;  DCS data tools
;
; AUTHOR:
;  Yiming Qiu
;  NIST Center for Neutron Research
;  100 Bureau Drive, Gaithersburg, MD 20899-6102
;  United States
;  yiming.qiu@nist.gov
;  July, 2007
;
; LICENSE:
;  The software in this file is written by an employee of
;  National Institute of Standards and Technology
;  as part of the DAVE software project.
;
;  The DAVE software package is not subject to copyright protection
;  and is in the public domain. It should be considered as an
;  experimental neutron scattering data reduction, visualization, and
;  analysis system. As such, the authors assume no responsibility
;  whatsoever for its use, and make no guarantees, expressed or
;  implied, about its quality, reliability, or any other
;  characteristic. The use of certain trade names or commercial
;  products does not imply any endorsement of a particular product,
;  nor does it imply that the named product is necessarily the best
;  product for the stated purpose. We would appreciate acknowledgment
;  if the DAVE software is used or if the code in this file is
;  included in another product.
;
;#######################################################################

pro dcs_cryst_align_e2lambda_event,event
    widget_control,/hourglass
    widget_control,event.handler,get_uvalue=state

    ;catch and ignore all errors in this program
    catch, myerror
    if myerror ne 0 then begin
       ok = dialog_message(dialog_parent=state.tlb,!error_state.msg,/error)
       catch,/cancel
       return
    end

    case event.id of
       state.eTxt:     dcs_cryst_align_e2lambda_update,state,/E
       state.lTxt:     dcs_cryst_align_e2lambda_update,state,/lambda
       state.plotBut:  launch_DCSToolsPerformance,event
       state.exitBut:  widget_control,event.handler,/destroy
       else:
    endcase
end

pro dcs_cryst_align_e2lambda_update,state,E=E,lambda=lambda
    eng = !values.f_nan & wav = !values.f_nan   ;initialization
    ;retrieve values from text boxes
    if keyword_set(E) then begin
       widget_control,state.eTxt,get_value=str
       tmp = dm_to_number(str)
       if finite(tmp) then begin
          eng = tmp
          if eng gt 0 then begin
             wav = sqrt(81.8042/eng)
             widget_control,state.lTxt,set_value=dm_to_string(wav)
          endif
       endif
    endif
    if keyword_set(lambda) then begin
       widget_control,state.lTxt,get_value=str
       tmp = dm_to_number(str)
       if finite(tmp) then begin
          wav = tmp
          if wav gt 0 then begin
             eng = 81.8042/wav/wav
             widget_control,state.eTxt,set_value=dm_to_string(eng)
          endif
        endif
    endif
    ;find out the energy resolution and flux
    if finite(eng) and finite(wav) then begin
       f_low  = interpol((*(state.flPtr))[1,*],(*(state.flPtr))[0,*],wav)/state.fmax*100.0
       f_med  = interpol((*(state.fmPtr))[1,*],(*(state.fmPtr))[0,*],wav)/state.fmax*100.0
       flux   = 'low:'+strtrim(string(0>(f_low)<100,format='(f6.2)'),2)+', '+$
                'med:'+strtrim(string(0>(f_med)<100,format='(f6.2)'),2)
       widget_control,state.fxTxt,set_value=flux
       er_low = 13220./wav^2.947
       er_med = 3948./wav^2.656
       eresl  = 'low:'+strtrim(string(er_low,format='(f9.1)'),2)+', '+$
                'med:'+strtrim(string(er_med,format='(f9.1)'),2)
       widget_control,state.erTxt,set_value=eresl
    endif else begin
       widget_control,state.fxTxt,set_value=''
       widget_control,state.erTxt,set_value=''
       if keyword_set(E) then      widget_control,state.lTxt,set_value=''
       if keyword_set(lambda) then widget_control,state.eTxt,set_value=''
    endelse
end

pro dcs_cryst_align_e2lambda_Exit,tlb
    widget_control,tlb,get_uvalue=state
    ptr_free,state.flPtr,state.fmPtr
    widget_control,tlb,/destroy
end

;E to Lambda conversion
pro dcs_cryst_align_e2lambda,event,lambda=lambda
    state={group_leader:    0L, $   ;group leader
           tlb:             0L, $   ;top level base
           eTxt:            0L, $   ;E text box
           lTxt:            0L, $   ;lambda text box
           erTxt:           0L, $   ;E resolution text box
           fxTxt:           0L, $   ;flux text box
           plotBut:         0L, $   ;start performance button
           exitBut:         0L, $   ;exit button
           fmax:            0e, $   ;maximum flux
           flPtr:    ptr_new(), $   ;low resl mode flux
           fmPtr:    ptr_new()  $   ;medium resl mode flux
          }
    registerName = 'dcs_cryst_align_e2lambda'
    if xregistered(registerName) then begin ;only allow one copy running
       id = LookupManagedWidget(registername)
       widget_control,id,iconify=0
       return
    endif
    if n_elements(lambda) eq 0 then lambda=!values.f_nan
    if n_elements(event) ne 0 then begin
       state.group_leader = event.top
       state.tlb = widget_base(title='Convert Energy to Wavelength',/col,kill_notify='dcs_cryst_align_e2lambda_Exit',$
                   tlb_frame_attr=1,group_leader=event.top)
    endif else $
       state.tlb = widget_base(title='Convert E to Wavelength',/col,kill_notify='dcs_cryst_align_e2lambda_Exit',$
                   tlb_frame_attr=1)
    row1        = widget_base(state.tlb,/row,/grid_layout)
    col1        = widget_base(row1,/col,/grid_layout)
    col2        = widget_base(row1,/col,/grid_layout)
    void        = widget_label(col1,value='E (meV)',/align_center)
    state.eTxt  = widget_text(col1,/editable,/all_events)
    void        = widget_label(col2,value='Lambda ('+string('c5'xb)+')',/align_center)
    state.lTxt  = widget_text(col2,/editable,/all_events,value=dm_to_string(lambda))
    void        = widget_label(col1,value='FWHM ('+string('b5'xb)+'eV)',/align_center)
    state.erTxt = widget_text(col1)
    void        = widget_label(col2,value='Flux (% of maximum)',/align_center)
    state.fxTxt = widget_text(col2)
    row2        = widget_base(state.tlb,/row,/align_center)
    state.plotBut = widget_button(row2,value='  Plot  ',sensitive=xregistered('dave',/noshow))
    void          = widget_label(row2,value= '   ')
    state.exitBut = widget_button(row2,value='  Exit  ')
    ;flux info
    state.flPtr = ptr_new([[1.25,357], [1.3,1799.4], [1.35,30630], [1.36,44460], [1.38,69780], [1.4,86595], [1.42,86820], [1.44,59370], [1.45,57570], [1.46,73200],$
                    [1.48,117900], [1.5,137498], [1.52,143280], [1.54,145380], [1.55,135750], [1.56,106530], [1.58,33600], [1.6,34440], [1.62,122010], [1.64,177630],$
                    [1.65,172500], [1.66,158280], [1.68,167760], [1.7,212850], [1.72,165270], [1.74,68310], [1.75,51450], [1.76,71100], [1.78,186870], [1.8,266385],$
                    [1.82,265980], [1.84,216420], [1.85,177360], [1.86,131010], [1.88,43860], [1.9,4959], [1.95,3192], [1.96,7602], [1.98,32730], [2,87915],$
                    [2.02,116220], [2.04,60900], [2.05,28993.5], [2.06,10188], [2.08,820.8], [2.1,167.22], [2.15,385.5], [2.2,196260], [2.25,541500], [2.3,720400],$
                    [2.4,803100], [2.5,875400], [2.6,905700], [2.7,912600], [2.8,902400], [2.9,898500], [3,901350], [3.1,888300], [3.2,876600], [3.25,855900],$
                    [3.27,826800], [3.3,760800], [3.33,390900], [3.35,346200], [3.38,731100], [3.4,831900], [3.42,839400], [3.45,846000], [3.5,861300], [3.6,870300],$
                    [3.7,877800], [3.8,859500], [3.85,844500], [3.9,832200], [3.95,807000], [3.98,781200], [4,763800], [4.03,734700], [4.05,787800], [4.08,850500],$
                    [4.1,872400], [4.15,857400], [4.2,847200], [4.3,828900], [4.4,798300], [4.5,765900], [4.6,738000], [4.62,719400], [4.65,716400], [4.67,746400],$
                    [4.7,780000], [4.72,777900], [4.75,772800], [4.78,766200], [4.8,763950], [5,721650], [5.5,645900], [6,558000], [6.3,508200], [6.4,487500],$
                    [6.45,466500], [6.5,453600], [6.55,405600], [6.6,297180], [6.65,33900], [6.7,24780], [6.75,279690], [6.8,384600], [6.85,402600], [6.9,413700],$
                    [6.95,406200], [7,404850], [7.5,349800], [8,292040], [9,199770], [10,142250], [12,69260], [14,34890], [15,25104], [16,17742],$
                    [18,9459], [20,5184]])
    state.fmPtr = ptr_new([[1.35,3894], [1.4,11083.5], [1.45,4354.5], [1.5,17070], [1.55,17355], [1.6,1409.55], [1.65,22560], [1.7,27450], [1.75,3330], [1.8,34275],$
                    [1.85,22425], [1.9,207.45], [1.95,322.2], [2,11304], [2.05,1801.5], [2.1,130.185], [2.15,150.9], [2.2,28005], [2.25,72900], [2.3,92160],$
                    [2.4,102165], [2.5,112905], [2.6,117765], [2.7,120540], [2.8,117645], [2.9,119475], [3,122535], [3.2,117090], [3.25,114390], [3.3,103335],$
                    [3.35,31950], [3.4,112305], [3.5,116445], [3.6,120810], [3.7,120810], [3.8,119025], [3.85,118320], [3.9,115035], [3.95,111540], [4,105735],$
                    [4.05,106050], [4.1,123960], [4.15,121905], [4.2,121185], [4.3,117795], [4.4,113085], [4.6,103245], [4.8,110460], [5,103703], [5.5,91095],$
                    [6,77842.5], [6.5,59820], [6.6,36780], [6.7,952.5], [6.8,49335], [6.9,53190], [7,54127.5], [7.5,43965], [8,36540], [9,24315],$
                    [10,16657.5], [12,7679.25], [14,3831], [16,1926], [18,1006.58], [20,554.55]])
    state.fmax = max((*(state.flPtr))[1,*])
    dm_center_kid,state.tlb,state.group_leader
    widget_control,state.tlb,/realize
    dcs_cryst_align_e2lambda_update,state,/lambda
    widget_control,state.tlb,set_uvalue=state
    xmanager,registerName,state.tlb,cleanup='dcs_cryst_align_e2lambda_Exit',/no_block
end

;plotting
pro dcs_cryst_align_plot,pState
    xtit = ['Rotation (!Eo!n)','Height (mm)','X-translation (mm)','Y-translation (mm)','X-tilt (!Eo!n)','Y-tilt (!Eo!n)',$
            'File Index','T (K)']
    xtit = xtit[(*pState).xaxis]
    ytit = 'Neutron Counts/'+strtrim(string((*pState).duration,format='(i5)'),2)+' Secs'
    tit  = (*pState).title  ;reduce the size
    x    = reform((*(*pState).xPtr)[(*pState).xaxis,*])
    mon  = reform((*(*pState).xPtr)[10,*]) &  mon_f = mon[0]
    n_file = n_elements((*(*pState).yPtr)[0,*])
    y = fltarr(n_file) & dy = y
    for i=0,n_file-1 do begin
        y[i]  = total((*(*pState).yPtr)[*(*pState).detPtr,i])
        dy[i] = sqrt(total((*(*pState).dyPtr)[*(*pState).detPtr,i]^2))
    endfor
    if ((*pState).xaxis eq 7) and (n_file gt 1) then begin  ;combine data with same T_set
       T_set = reform((*(*pState).xPtr)[9,*])
       ind   = sort(T_set)
       ndat  = n_elements(T_set)
       x     = x[ind]
       y     = y[ind]
       mon   = mon[ind]
       T_set = T_set[ind]
       y     = y*mon/mon_f
       ind   = uniq(T_set)
       nun   = n_elements(ind)
       dy    = fltarr(nun)
       j     = 0L
       for i=0L,nun-1 do begin
           xtmp = 0.
           ytmp = 0.
           mtmp = 0.
           itmp = 0L
           while(T_set[j] eq T_set[ind[i]]) do begin
                xtmp = xtmp+x[j]
                ytmp = ytmp+y[j]
                mtmp = mtmp+mon[j]
                itmp = itmp+1L
                j    = j+1L
                if j eq ndat then break
           endwhile
           x[i]  = xtmp/(itmp>1)
           y[i]  = ytmp/mtmp*mon_f
           dy[i] = sqrt(ytmp)/mtmp*mon_f
       endfor
       x = x[0:nun-1]
       y = y[0:nun-1]
    endif
    xmax = max(x,min=xmin)
    ymax = max(y+dy,imax) & imax = imax[0]
    ymin = min(y-dy)
    xran = [xmin,xmax]+[-0.05,0.05]*(xmax-xmin+abs(xmax+1)/10.)
    yran = [ymin,ymax]+[-0.05,0.05]*(ymax-ymin)
    (*pState).plotWin->setproperty,xdat=x,ydat=y,yerr=dy,xran=xran,yran=yran,$
       title=tit,xtit=xtit,ytit=ytit,psym='circle',linestyle='no line',/hidelegend,legend='data'
    (*pState).plotWin->remove_text
    index = uniq(x,sort(x))
    if ((*pState).fit eq 1) and (n_elements(index) ne 1) and (n_elements(y) gt 4)then begin
       fwhm = total(y-ymin)/(ymax-ymin)*abs(x[1]-x[0])
       estimates    = [ymax,x[imax],fwhm/sqrt(8.*alog(2)),ymin]
       estimates[2] = estimates[2]/2.
       result = gaussfit(x,y,params,estimates=estimates,nterms=n_elements(estimates))
       xx = xmin+findgen(150)/149*(xmax-xmin+2)-1
       yy = params[0]*exp(-0.5*((xx-params[1])/params[2])^2)+params[3]
       (*pState).plotWin->add_plot,xx,yy,linestyle='solid',color='red',legend='gauss fit'
       (*pState).plotWin->add_text,'peak!l !n='+strtrim(string(params[1],format='(f8.2)'),2)+'!C'+$
                'fwhm='+strtrim(string(params[2]*sqrt(8.*alog(2)),format='(f9.3)'),2),0.04,0.85,$
                fontsize=10
    endif
    (*pState).plotWin->draw
end

;calculate the angle between hkl1 and hkl2
function dcs_cryst_align_angle,hkl1,hkl2,dir_latt_parm
    abc_r = dm_basis_r(dir_latt_parm)
    q1 = abc_r#hkl1
    q2 = abc_r#hkl2
    return,acos(q1##transpose(q2)/norm(q1)/norm(q2))/!dtor
end

;calculate q length from direct crystal lattice parameters
function dcs_cryst_align_qlength,hkl,cryst
    ddtor   = !dpi/(180d)
    sinabg  = sin(cryst[3:5]*ddtor)
    cosabg  = cos(cryst[3:5]*ddtor)
    v       = cryst[0]*cryst[1]*cryst[2]*sqrt((1d)+(2d)*cosabg[0]*cosabg[1]*cosabg[2]-total(cosabg[0:2]^2))
    rec_abc = (2d)*!dpi/v*[cryst[1]*cryst[2]*sinabg[0],cryst[0]*cryst[2]*sinabg[1],cryst[0]*cryst[1]*sinabg[2]]
    tmp     = (4d)*!dpi*!dpi/v/v*cryst[0]*cryst[1]*cryst[2]
    aadotbb = tmp*cryst[2]*(cosabg[0]*cosabg[1]-cosabg[2])
    bbdotcc = tmp*cryst[0]*(cosabg[1]*cosabg[2]-cosabg[0])
    ccdotaa = tmp*cryst[1]*(cosabg[2]*cosabg[0]-cosabg[1])
    return, sqrt(total((hkl*rec_abc)^2)+(2d)*(hkl[0]*hkl[1]*aadotbb+hkl[1]*hkl[2]*bbdotcc+hkl[2]*hkl[0]*ccdotaa))
end

pro dcs_cryst_align_event,event
    compile_opt IDL2    ;,strictarrsubs
                        ;idl2=defint32,strictarr
                        ;strictarr:   [] to index array
                        ;strictarrsubs: error when out-of-range indices,IDL5.6 or after
    WIDGET_CONTROL,/HOURGLASS
    widget_control,event.handler,get_uvalue=pState

    if (*pState).enotebook then $
       IF DAVE_SET_FOCUS(EVENT) THEN RETURN

    widget_type = widget_info(event.id,/type)
    case widget_type of
       8:    tmp_value = widget_info(event.id,/droplist_select)  ;droplist
       12:   tmp_value = event.index                             ;combobox
       else: tmp_value = 0
    endcase

    ;catch and ignore all errors in this program
    catch, myerror
    if myerror ne 0 then begin
       ok = dialog_message(dialog_parent=(*pState).tlb,!error_state.msg,/error)
       catch,/cancel
       return
    end

    case (event.id) of
       (*pState).tlb:    widget_control,event.id,scr_xsize=(*pState).geom[0],scr_ysize=(*pState).geom[1]
       (*pState).loadBut: begin
          (*pState).fit = 0
          (*pState).filesel->getproperty,path=open_path,file=open_file
          if n_elements(open_file) eq 0 then begin
             ok = dialog_message('Please select DCS files first.',/error,dialog_parent=event.handler)
             break
          endif
          e_range = (*pState).e_ran
          dm_load_dcsdata,open_path,open_file,ang=x,qty=y,dqty=dy,t_chan=(*pState).t_chan,e_range=e_range,$
              parent=(*pState).tlb,duration=duration,/raw,/eadjust,montype=0
          if n_elements(e_range) ne 0 then (*pState).e_ran = e_range
          *(*pState).xPtr    = temporary(x)
          *(*pState).yPtr    = temporary(y)
          *(*pState).dyPtr   = temporary(dy)
          (*pState).duration = fix(duration,type=3)
          dcs_cryst_align_plot,pState
         end
       (*pState).fitBut:  begin
          if n_elements(*(*pState).yPtr) le 4 then break
          (*pState).fit = 1
          dcs_cryst_align_plot,pState
         end
       (*pState).saveBut: begin
          if n_elements(*(*pState).xPtr) eq 0 then break
          (*pState).plotWin->saveas,'ps'
         end
       (*pState).dataBut:  begin
          if n_elements(*(*pState).xPtr) eq 0 then break
          file = dm_choose_file('txt',/write,dialog_parent=(*pState).tlb,title='Save plot data as an ASCII file.')
          if strlen(file) eq 0 then break
          openw,unit,file,/get_lun,error=openerr
          if openerr ne 0 then begin
             ok = dialog_message("Can't write in "+file[0],/error,dialog_parent=(*pState).tlb)
          endif else begin
             n_file = n_elements((*(*pState).yPtr)[0,*])
             widget_control,/hourglass  ;busy signal
             x   = reform((*(*pState).xPtr)[(*pState).xaxis,*])
             mon = reform((*(*pState).xPtr)[10,*]) &  mon_f = mon[0]
             n_file = n_elements((*(*pState).yPtr)[0,*])
             y = fltarr(n_file) & dy = y
             for i=0L,n_file-1L do begin
                 y[i]  = total((*(*pState).yPtr)[*(*pState).detPtr,i])
                 dy[i] = sqrt(total((*(*pState).dyPtr)[*(*pState).detPtr,i]^2))
             endfor
             if ((*pState).xaxis eq 7) and (n_file gt 1) then begin  ;combine data with same T_set
                T_set = reform((*(*pState).xPtr)[9,*])
                ind   = sort(T_set)
                ndat  = n_elements(T_set)
                x     = x[ind]
                y     = y[ind]
                mon   = mon[ind]
                T_set = T_set[ind]
                y     = y*mon/mon_f
                ind   = uniq(T_set)
                nun   = n_elements(ind)
                dy    = fltarr(nun)
                j     = 0L
                for i=0L,nun-1 do begin
                    xtmp = 0.
                    ytmp = 0.
                    mtmp = 0.
                    itmp = 0L
                    while(T_set[j] eq T_set[ind[i]]) do begin
                        xtmp = xtmp+x[j]
                        ytmp = ytmp+y[j]
                        mtmp = mtmp+mon[j]
                        itmp = itmp+1L
                        j    = j+1L
                        if j eq ndat then break
                    endwhile
                    x[i]  = xtmp/(itmp>1)
                    y[i]  = ytmp/mtmp*mon_f
                    dy[i] = sqrt(ytmp)/mtmp*mon_f
                endfor
                x = x[0:nun-1]
                y = y[0:nun-1]
             endif
             for i=0L,n_elements(x)-1L do $
                 printf,unit,x[i],' ',y[i],' ',dy[i]
             free_lun,unit
          endelse
         end
       (*pState).printBut: begin
          if n_elements(*(*pState).xPtr) eq 0 then break
          (*pState).plotWin->saveas,'printer'
         end
       (*pState).exitBut: begin
          widget_control,event.handler,/destroy
          return
         end
       (*pState).detList: begin
          old = *(*pState).detPtr
          case tmp_value of
              0: begin
                 *(*pState).detPtr = indgen((*pState).ndet)
                 (*pState).det_tit = 'det=all'
                 end
              1: begin
                 n_index = n_elements(*(*pState).detPtr)
                 if n_index eq (*pState).ndet then default='1:'+strtrim(string((*pState).ndet),2) $
                 else default=string(*(*pState).detPtr+1,format='('+string(n_index)+'(i3,","))')
                 if strmid(default,0,/reverse_offset) eq ',' then default=strmid(default,0,strlen(default)-1)
                 tmp = dm_dialog_input(dialog_parent=event.handler,'Det num:',default=default,title='Specifying detectors',$
                       info='Detector number starts from 1. eg:371,374 or 370:374 c',xsize=200)
                 (*pState).det_tit = 'det='+tmp
                 if strlen(tmp) ne 0 then begin
                    tmpstrs = strsplit(tmp,',',/extract,count=count1)
                    for ii=0,count1-1 do begin
                        status = ''
                        index  = strpos(tmpstrs[ii],'c',/reverse_search)
                        if index[0] ne -1 then status='c'
                        index  = strpos(tmpstrs[ii],'C',/reverse_search)
                        if index[0] ne -1 then status='c'
                        index  = strpos(tmpstrs[ii],'u',/reverse_search)
                        if index[0] ne -1 then status='u'
                        index  = strpos(tmpstrs[ii],'U',/reverse_search)
                        if index[0] ne -1 then status='u'
                        index  = strpos(tmpstrs[ii],'l',/reverse_search)
                        if index[0] ne -1 then status='l'
                        index  = strpos(tmpstrs[ii],'L',/reverse_search)
                        if index[0] ne -1 then status='l'

                        ind0  = strsplit(tmpstrs[ii],'[a-zA-Z '+string(9b)+',:-]',/regex)
                        ind1  = strsplit(tmpstrs[ii],':-')
                        index = -1L
                        for i=1,n_elements(ind1)-1 do begin
                            ind2 = where(ind0 eq ind1[i],count)
                            if count ne 0 then index=[index,ind2]
                        endfor
                        n_index = n_elements(index)
                        if n_index ne 1 then index=index[1:n_index-1]-1
                        tmp = strsplit(tmpstrs[ii],'[a-zA-Z '+string(9b)+',:-]',/regex,/extract)
                        det_tmp = -1L
                        for i=0,n_elements(tmp)-1 do begin
                            tmp1 = dm_to_number(tmp[i],/float)
                            if finite(tmp1) then begin
                               ind0 = where(index eq i,count)
                               if count ne 0 then begin
                                  i = i+1
                                  tmp2 = dm_to_number(tmp[i],/float)
                                  if finite(tmp2) then det_tmp=[det_tmp,tmp1+findgen(abs(tmp2-tmp1)+1)] $
                                  else det_tmp=[det_tmp,tmp1]
                               endif else det_tmp=[det_tmp,tmp1]
                            endif
                        endfor
                        index = where(det_tmp gt 0,count)
                        if count ne 0 then begin
                           det_tmp=fix(det_tmp[index],type=3)-1
                           case status of
                                'c': det_tmp = dm_common(det_tmp,dm_mask_detbanks(0),/no_copy)
                                'l': det_tmp = dm_common(det_tmp,dm_mask_detbanks(-1),/no_copy)
                                'u': det_tmp = dm_common(det_tmp,dm_mask_detbanks(1),/no_copy)
                                else:det_tmp = dm_common(det_tmp,indgen((*pState).ndet),/no_copy)
                           endcase
                           if det_tmp[0] ne -1 then begin
                              if n_elements(dets) eq 0 then $
                                 dets = det_tmp $
                              else $
                                 dets = [dets,det_tmp]
                           endif
                        endif
                    endfor
                 endif
                 if (n_elements(dets) eq 0) or (n_elements(dets) eq (*pState).ndet) then begin
                    *(*pState).detPtr = indgen((*pState).ndet)
                    (*pState).det_tit = 'det=all'
                    dm_set_droplist,(*pState).detList,select=0
                 endif else begin
                    *(*pState).detPtr = dets
                    tth = strtrim(string(mean(abs((*(*pState).tthPtr)[dets])),format='(f9.2)'),2)
                    if n_elements(dets) gt 1 then tth='(<2th>='+tth+'!Eo!n)' $
                    else tth='(2th='+tth+'!Eo!n)'
                    (*pState).det_tit = (*pState).det_tit+tth
                 endelse
                 end
              2: begin
                 *(*pState).detPtr = dm_mask_detbanks(0)
                 (*pState).det_tit = 'det=central'
                 end
              3: begin
                 *(*pState).detPtr = dm_mask_detbanks(1)
                 (*pState).det_tit = 'det=upper'
                 end
              4: begin
                 *(*pState).detPtr = dm_mask_detbanks(-1)
                 (*pState).det_tit = 'det=lower'
                 end
              else:begin
                 *(*pState).detPtr = tmp_value-5
                 tth = strtrim(dm_to_string(abs((*(*pState).tthPtr)[tmp_value-5])),2)
                 (*pState).det_tit = 'det='+strtrim(string(tmp_value-4),2)+'(2th='+tth+'!Eo!n)'
                 end
          endcase
          (*pState).title = (*pState).det_tit+', '+(*pState).tch_tit
          n_old = n_elements(old)
          n_new = n_elements(*(*pState).detPtr)
          comm  = dm_common(old,*(*pState).detPtr)
          same  = (n_old eq n_new) and (n_old eq n_elements(comm)) and (comm[0] ne -1)
          old   = -1 & comm = -1
          if (n_elements(*(*pState).yPtr) ne 0) and (1-same) then begin   ;reload data
             dcs_cryst_align_plot,pState
          endif
         end
       (*pState).timeList:   if ((*pState).t_chan ne tmp_value) or (tmp_value eq 2) then begin
          (*pState).t_chan = tmp_value
          if (*pState).t_chan eq 0 then (*pState).tch_tit='tchan=all' $
          else if (*pState).t_chan eq 1 then (*pState).tch_tit='tchan=elastic' $
          else begin
             default = dm_to_string((*pState).e_ran)
             new     = dm_dialog_input(['from:','to:'],title='Energy range',xsize=100,default=default,/float,$
                       dialog_parent=(*pState).tlb,info='In units of meV')
             (*pState).t_chan = 1
             (*pState).e_ran[*] = !values.f_nan
             if total(finite(new)) eq 2 then $
                if new[1] gt new[0] then begin
                   (*pState).e_ran   = new
                   (*pState).t_chan  = 2
                   (*pState).tch_tit = 'E=['+dm_to_string(new[0])+','+dm_to_string(new[1])+']meV'
                endif
             if (*pState).t_chan eq 1 then begin
                ok = dialog_message('Invalid energy range. Elastic peak option is used.',dialog_parent=(*pState).tlb)
                dm_set_droplist,event.id,select=1
                (*pState).tch_tit = 'tchan=elastic'
             endif
          endelse
          (*pState).title = (*pState).det_tit+', '+(*pState).tch_tit
          (*pState).filesel->getproperty,file=open_file
          if n_elements(open_file) ne 0 then begin ;reload data
             tlb  = (*pState).tlb
             load = (*pState).loadBut
             dcs_cryst_align_event,{WIDGET_BUTTON,ID:load,TOP:tlb,HANDLER:tlb,SELECT:1}
             return
          endif
         endif
       (*pState).xaxisList:begin
          (*pState).xaxis = tmp_value
          (*pState).fit   = 0
          if n_elements(*(*pState).xPtr) ne 0 then dcs_cryst_align_plot,pState
         end
       (*pState).calctthBut:begin
          default = (float(*(*pState).infoPtr))[0:9]
          info    = dm_dialog_input(['lambda','a','b','c','alpha','beta','gamma','H','K','L'],$
                    default=default,/float,dialog_parent=(*pState).tlb,xsize=80,title='Calculate 2theta')
          if total(finite(info)) eq n_elements(default) then begin
             (*(*pState).infoPtr)[0:9] = double(info)
             ki = 2.*!dpi/info[0]
             q  = dcs_cryst_align_qlength(info[7:9],info[1:6])
             if q/2. lt ki then begin
                mesg = ''
                tth  = asin(q/2./ki)*2/!dtor
                idet = where(*(*pState).tthPtr ge tth,count)
                if count ne 0 then begin
                   if idet[0]-1 gt 0 then begin
                      if abs((*(*pState).tthPtr)[idet[0]-1]-tth) lt abs((*(*pState).tthPtr)[idet[0]]-tth) then $
                         idet[0] = idet[0]-1
                   endif
                   tmp  = [idet[0],(912<(idet[0]+1)),(0>(idet[0]-1))]
                   idet = dm_common(tmp,dm_mask_detbanks(0),count)
                   if count eq 0 then $
                      mesg = 'The corresponding detector number is '+dm_to_string(tmp[0]+1,/int)+'.'  $
                   else $
                      mesg = 'The corresponding central bank detector number is '+dm_to_string(idet[0]+1,/int)+'.'
                endif else $
                   mesg = 'It is beyond the detector range.'
                tth  = strtrim(string(tth,format='(f9.2)'),2)+string('b0'xb)
                mesg = ['Two theta angle for ['+dm_to_string(info[7])+','+dm_to_string(info[8])+','+$
                       dm_to_string(info[9])+'] is '+tth+'.',mesg]
             endif else begin
                lmax = strtrim(string(4.*!dpi/q,format='(f9.2)'),2)+string('c5'XB)
                mesg = 'Incident wavelength must be smaller than '+lmax+' to view ['+dm_to_string(info[7])+','+$
                        dm_to_string(info[8])+','+dm_to_string(info[9])+'] reflection.'
             endelse
             ok = dialog_message(mesg,/info,dialog_parent=(*pState).tlb,title='Result')
          endif
         end
       (*pState).calclattBut:begin
          default = strpos((*pState).det_tit,'!9%')
          if default[0] ne -1 then begin
             default = strsplit((*pState).det_tit,'=',/extract)
             if n_elements(default) ge 3 then $
                default = dm_to_number(default[2],/float)    $
             else default = !values.f_nan
          endif else $
             default = !values.f_nan
          default = [default,float((*(*pState).infoPtr)[[0,4,5,6,7,8,9]])]
          info = dm_dialog_input(['2theta','lambda','alpha','beta','gamma','H','K','L'],title='calculate',$
                 info='assume a=b=c',default=default,/float,dialog_parent=(*pState).tlb,xsize=80)
          if total(finite(info)) eq 8 then begin
             info = double(info)
             qlen = dcs_cryst_align_qlength(info[5:7],[1.0,1.0,1.0,info[2:4]])
             a    = qlen*info[1]/4./!dpi/sin(info[0]/2*!dtor)
             ok   = dialog_message("The lattice paramter is "+dm_to_string(float(a))+string('c5'xB)+' assuming a=b=c.',$
                  /info,dialog_parent=(*pState).tlb,title='Result')
          endif
         end
       (*pState).calca2But:begin
          default = strpos((*pState).det_tit,'!9%')
          if default[0] ne -1 then begin
             default = strsplit((*pState).det_tit,'=',/extract)
             if n_elements(default) ge 3 then $
                default = dm_to_number(default[2],/float)    $
             else default = !values.f_nan
          endif else $
             default = !values.f_nan
          tth = dm_dialog_input('two theta',default=default,info='A2=90-theta',xsize=80,$
                dialog_parent=(*pState).tlb,title='Calculate A2')
          if finite(tth) then begin
             mesg = 'Set A2 to '+strtrim(string(90.-tth/2.,format='(f9.2)'),2)+string('b0'xb)+'.'
             ok = dialog_message(mesg,/info,dialog_parent=(*pState).tlb,title='Result')
          endif
         end
       (*pState).calca2Rot:begin
          default = (float(*(*pState).infoPtr))[0:12]
          info    = dm_dialog_input(['lambda','a','b','c','alpha','beta','gamma','H1','K1','L1','H2','K2','L2'],$
                    default=default,/float,dialog_parent=(*pState).tlb,xsize=80,title='Calculate A2 Rotation',$
                    info='1->2 counterclockwise')
          if total(finite(info)) eq n_elements(default) then begin
             (*(*pState).infoPtr)[0:12] = double(info)
             ki = 2.*!dpi/info[0]
             q1 = dcs_cryst_align_qlength(info[7:9],info[1:6])
             q2 = dcs_cryst_align_qlength(info[10:12],info[1:6])
             if max([q1,q2])/2. lt ki then begin
                mesg = ''
                tth1 = asin(q1/2./ki)*2/!dtor
                tth2 = asin(q2/2./ki)*2/!dtor
                idet = where(*(*pState).tthPtr ge max([tth1,tth2]),count)
                if count ne 0 then begin
                   abet = dcs_cryst_align_angle(info[7:9],info[10:12],info[1:6]) ;angle between q1 and q2
                   ang1 = 90.0-tth1/2.
                   ang2 = 90.0-tth2/2.-abet
                   mesg = 'A2 needs to rotate '+strtrim(string(ang2-ang1,format='(f9.2)'),2)+string('b0'xb)+'.'
                endif else $
                   mesg = 'The vectors are beyond the detector range.'
             endif else begin
                lmax = strtrim(string(4.*!dpi/(max([q1,q2])),format='(f9.2)'),2)+string('c5'XB)
                mesg = 'Incident wavelength must be smaller than '+lmax+'.'
             endelse
             ok = dialog_message(mesg,/info,dialog_parent=(*pState).tlb,title='Result')
          endif
         end
       (*pState).e2lambdaBut: dcs_cryst_align_e2lambda,event,lambda=float((*(*pState).infoPtr)[0])
       (*pState).tools: widget_displaycontextmenu,event.id,event.x,event.y,(*pState).contextMenu
       else:
    endcase
end

pro dcs_cryst_align_Exit,tlb
    widget_control,tlb,get_uvalue=pState,/no_copy
    obj_destroy,(*pState).filesel
    obj_destroy,(*pState).plotWin
    cd,(*pState).root
    ptr_free,(*pState).detPtr,(*pState).xPtr,(*pState).yPtr,(*pState).dyPtr,(*pState).tthPtr,(*pState).infoPtr
    ptr_free,pState
    widget_control,tlb,/destroy
end

;initparms=[a,b,c,alpha,beta,gamma,uh,uk,ul]
pro dcs_cryst_align,event,initparms
    state={group_leader:                 0L, $  ;group leader
           tlb:                          0L, $  ;top level base
           geom:                    [0L,0L], $  ;geometry
           filesel:               obj_new(), $  ;file selector
           root:                         '', $  ;starting directory
           infoPtr: ptr_new(/allocate_heap), $  ;store the info for calculating 2theta
                                             $  ;[lambda,a,b,c,alpha,beta,gamma,h1,k1,l1,h2,k2,l2]
           contextMenu:                  0L, $  ;context sensitive menu base
           calctthBut:                   0L, $  ;calculate 2theta button
           calclattBut:                  0L, $  ;calculate cubic lattice parameter
           calca2But:                    0L, $  ;calculate A2 button
           calca2Rot:                    0L, $  ;calculate A2 rotation button
           e2lambdaBut:                  0L, $  ;E to lambda conversion button
           loadBut:                      0L, $  ;load file button
           fitBut:                       0L, $  ;gauss fit button
           fit:                          0b, $  ;flag for fit
           saveBut:                      0L, $  ;save plot button
           dataBut:                      0L, $  ;save data button
           printBut:                     0L, $  ;print plot button
           exitBut:                      0L, $  ;exit button
           ndet:                         0L, $  ;number of detectors
           detList:                      0L, $  ;detector droplist
           tools:                        0L, $  ;base for the tools context sensitive menu
           tthPtr:  ptr_new(/allocate_heap), $  ;save the two theta value
           detPtr:  ptr_new(/allocate_heap), $  ;pointer to detector array
           timeList:                     0L, $  ;time channel droplist
           t_chan:                       1s, $  ;0-all 1-elastic 2-specify
           e_ran:                   [0e,0e], $  ;energy range for t_chan=2
           xaxisList:                    0L, $  ;x-axis droplist
           xaxis:                        0s, $  ;xaxis number
           xPtr:    ptr_new(/allocate_heap), $  ;pointer to x data
           yPtr:    ptr_new(/allocate_heap), $  ;pointer to y data
           dyPtr:   ptr_new(/allocate_heap), $  ;pointer to dy data
           plotWin:               obj_new(), $  ;draw widget
           duration:                     0L, $  ;duration of measurement in secs
           det_tit:                      '', $  ;output detector string
           tch_tit:                      '', $  ;output time channel string
           title:                        '', $  ;plot title
           workDir:                      '', $  ;work directory
           datDir:                       '', $  ;data directory
           enotebook:                    0b  $  ;use Rob's enotebook
    }
    state.e_ran  = [!values.f_nan,!values.f_nan]
    registerName = 'dcs_cryst_align'
;   if xregistered(registerName) then return  ;only allow one copy to be running at one time
    cd,current = current
    state.root = current
    state.workDir = current
    state.datDir  = current

    *state.infoPtr= dblarr(13) & (*state.infoPtr)[*] = !values.d_nan
    if (n_elements(initparms) le 12) and (n_elements(initparms) ne 0) then $
       (*state.infoPtr)[1:n_elements(initparms)] = initParms
    ;get the detector information
    alldetPtr     = dm_load_detpos(/dcs)
    *state.tthPtr = (*alldetPtr).two_theta
    if n_elements(event) ne 0 then begin
       state.group_leader = event.top
       state.tlb  = widget_base(title='DCS Single Crystal Alignment Utility',/row,kill_notify='dcs_cryst_align_Exit',/tlb_size_event,$
         group_leader=event.top,map=0)
    endif else $
       state.tlb  = widget_base(title='DCS Single Crystal Alignment Utility',/row,kill_notify='dcs_cryst_align_Exit',/tlb_size_event,$
         map=0)
    left = widget_base(state.tlb,/col)
    state.filesel = obj_new('dm_filesel',parent=left,xsize=160,ysize=20,/frame,filter='.dcs')
    defsysv,'!dave_defaults',exists=exists
    if exists then $
       if ptr_valid(!dave_defaults) then begin
         state.datDir  = (*!dave_defaults).datDir
         state.workDir = (*!dave_defaults).workDir
       endif
    state.filesel->set_path,state.datDir
    if !version.os_family eq 'unix' then $
       center = widget_base(state.tlb,/grid_layout,/col,/align_bottom,xpad=0,ypad=0,space=0) $
    else $
       center = widget_base(state.tlb,/grid_layout,/col,/align_bottom)

    ;construct the detector droplist values
    state.ndet = n_elements((*alldetPtr).two_theta)
    tmp   = string((*alldetPtr).two_theta,format='(f9.2)')
    tmp   = string(indgen(state.ndet)+1,format='(i3)')+tmp
    index = dm_mask_detbanks(0)
    tmp[index] = tmp[index]+'  c' & index=-1  ;central bank
    index = dm_mask_detbanks(1)
    tmp[index] = tmp[index]+'  u' & index=-1  ;upper bank
    index = dm_mask_detbanks(-1)
    tmp[index] = tmp[index]+'  l' & index=-1  ;lower bank
    tmp   = ['All','Specify','Central Bank','Upper Bank','Lower Bank',tmp]
    ptr_free,alldetPtr
    *state.detPtr  = indgen(state.ndet)  ;default all detectors
    state.det_tit  = 'det=all'
    state.tch_tit  = 'tchan=elastic'
    state.title    = state.det_tit+', '+state.tch_tit
    void = widget_label(center,value='Det. Number')
    state.detList  = dm_widget_droplist(center,value=tmp)
    void = widget_label(center,value='Time Channel')
    state.timeList = dm_widget_droplist(center,value=['All','Elastic','Specify'],select=1)
    void = widget_label(center,value='X-axis')
    state.xaxisList= dm_widget_droplist(center,value=['Rotation','Height','X-translation','Y-translation','X-tilt','Y-tilt','File Index','Temperature'])
    state.tools    = widget_base(center,/row,/context_events,/align_center)
    void    = widget_label(state.tools,value='Tools',sensitive=0)

    ;context sensitive tools menu
    state.contextMenu = widget_base(state.tools,/context_menu)
    state.calctthBut  = widget_button(state.contextMenu,value='Calculate 2theta')
    state.calclattBut = widget_button(state.contextMenu,value='Calculate Lattice Parameter')
    state.calca2But   = widget_button(state.contextMenu,value='Calculate A2 Reset Value')
    state.calca2Rot   = widget_button(state.contextMenu,value='Calculate A2 Rotation Angle')
    state.e2lambdaBut = widget_button(state.contextMenu,value='E to Wavelength Conversion',/separator)

    state.loadBut  = widget_button(center,value='Load Data')
    state.fitBut   = widget_button(center,value='Gauss Fit')
    state.saveBut  = widget_button(center,value='Save Plot')
    state.dataBut  = widget_button(center,value='Save Data')
    state.printBut = widget_button(center,value='Print')
    state.exitBut  = widget_button(center,value='Exit')

    geom  = widget_info(state.exitBut,/geometry)
    widget_control,state.detList,xsize=geom.scr_xsize
    widget_control,state.timeList,xsize=geom.scr_xsize
    widget_control,state.xaxisList,xsize=geom.scr_xsize

    geom  = widget_info(left,/geometry)
    ysize = geom.scr_ysize
    xsize = fix(1.3*ysize,type=3)
    state.plotWin  = obj_new('dm_plot',xsize=xsize,ysize=ysize,widgetbase=state.tlb,group_leader=state.tlb,$
        /compound,background='black')

    dm_center_kid, state.tlb   ;centering the top level base
    widget_control,state.tlb,/realize,/map

    geom  = widget_info(state.tlb,/geometry)
    state.geom = [geom.scr_xsize,geom.scr_ysize]

    pState = ptr_new(state,/no_copy)
    widget_control,(*pState).tlb,set_uvalue=pState

    if xregistered('dave',/noshow) and (n_elements(event) ne 0) then begin
       (*pState).enotebook = 1b
       RET = DAVE_SET_FOCUS((*pState).tlb)
    endif

    xmanager,registerName,(*pState).tlb,cleanup='dcs_cryst_align_Exit',/no_block
end