; Time-stamp: <2000-11-15 14:30:07 dph> ; MIT Directory: ~dph/h3/CXC/ARD/pro ; File: plt_osip.pro ; Author: D. Huenemoerder ; Original version: 2000.07.20 ; revised, 2000.08.20 ;==================================================================== ; visualize the OSIP (Order Sorting, Integrated Probability) table ; Given an event list's chip, chipx, chipy, tg_mlam, and an order, ; read the osip table, compute the order limits using empirical ; mapping of tg_mlam to chipx, chipy, and chip. ; Return a sorted array of x,y_lo, y_hi, where x is tg_mlam, and y are ; the limits (in units of real-valued diffraction order) ; Input arrays must be pre-selected for grating (e.g., heg or meg, via ; tg_part), and tg_srcid, if multiple sources. FUNCTION vis_osip, t_osip, chip_id, chipx, chipy, tg_mlam, order, fudge np = n_params() IF np LT 5 THEN BEGIN message, /inf, 'Usage: ' print, 'xy_struct = vis_osip(f_osip, chip_id, chipx, chipy, tg_mlam[, order[,fudge]])' print, ' t_osip = structure, OSIP ARD, as returned by mrdfits().' print, ' chip_id = array of chip_id (ccd_id) column from event list' print, ' chipx = array of chipx column from event list' print, ' chipy = array of chipy column from event list' print, ' tg_mlam = array of tg_mlam from event list' print, ' (chip_id, chipx, chipy, tg_mlam must be same length)' print, ' order = integer; order desired for limits; (default=1)' print, ' fudge = array of fudge-factors, one for each -S ccd' print, ' (Note: if fudge!=1, then FRACRESP is indef)' print, '' print, ' result is structure w/ tags: m_lo, m_hi; ' print, ' output array is commensurate w/ input arrays.' return, -1 ENDIF IF np EQ 5 THEN BEGIN order = 1 ; default to first order. fudge = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0] ; don't adjust limits ENDIF IF np EQ 6 THEN BEGIN fudge = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0] ; don't adjust limits ENDIF chip_id_S0 = 4 n_ev = n_elements(chip_id) dat = {m_lo:fltarr(n_ev), $ m_hi:fltarr(n_ev), $ fracresp:fltarr(n_ev)} FOR iev = 0L, n_ev-1 DO BEGIN c = chip_id[iev] cx = chipx[iev] cy = chipy[iev] mlam = tg_mlam[iev] fudge_factor = fudge[c-chip_id_S0] ; find entry in osip table; should ; result in one record. ll_osip = $ where($ (c EQ t_osip.ccd_id) AND $ (cx GE t_osip.chipx_min) AND $ (cx LE t_osip.chipx_max) AND $ (cy GE t_osip.chipy_min) AND $ (cy LE t_osip.chipy_max), $ n_osip) IF n_osip GT 1 THEN BEGIN message, /inf, 'Fatal error: non-unique osip-mapping.' print, 'chip,chipx,chipy=', c, cx, cy return, -1 ENDIF IF n_osip EQ 1 THEN BEGIN t = t_osip[ll_osip] ; interpolate table e_wish = !hc/abs(mlam) * order * 1000. e_lo = interpol_sort(t.energ_lo, t.energy, e_wish, /sil) e_hi = interpol_sort(t.energ_hi, t.energy, e_wish, /sil) dat.fracresp[iev] = interpol_sort(t.fracresp, t.energy, e_wish, /sil) ; dat.fracresp[iev] = interpol(t.fracresp, t.energy, e_wish) ; compute order limits dat.m_lo[iev] = abs(mlam)/(!hc/e_lo)/ 1000. / fudge_factor dat.m_hi[iev] = abs(mlam)/(!hc/e_hi)/ 1000. * fudge_factor ENDIF ENDFOR return, dat END ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO plt_osip, fevt, fosip, LBL=lbl, NEVT=n_evt, ORDERS=ords, FUDGE = fudge np = n_params() IF n_elements(lbl) EQ 0 THEN lbl = '' IF n_elements(n_evt) EQ 0 THEN n_evt = 0 ; means get them all IF n_elements(ords) EQ 0 THEN ords = [1] ; force to array IF n_elements(fudge) EQ 0 THEN fudge = replicate(1., 6) n_ords = n_elements(ords) ; set up output ; feroot = (str_sep( (reverse(str_sep(fevt,'/')))[0],'.fits'))[0] foroot = (str_sep( (reverse(str_sep(fosip,'/')))[0],'.fits'))[0] froot = feroot+lbl ;+'_'+foroot foo = ps_open() IF !ioflags.do_ps EQ 1 THEN BEGIN device, /port, /inch, /col, xsi = 6, ysi = 6, file = froot+'.ps' ENDIF ; read files: ; e=mrdfits(fevt,1,he, range = n_evt) t=mrdfits(fosip,1,ht) ; select good 1st order events: ; lgood = (e.status EQ 0) AND (e.grade LT 7) AND (e.grade NE 1) AND (e.grade NE 5) lm=where( (e.tg_part EQ 2) AND lgood ) lh=where( (e.tg_part EQ 1) AND lgood) order = 1 lm1=where( (abs(e.tg_m) EQ order) AND (e.tg_part EQ 2) AND lgood) lh1=where( (abs(e.tg_m) EQ order) AND (e.tg_part EQ 1) AND lgood) ; get osip limits per event for MEG (this is time-consuming): ; message, /inf, 'Computing MEG 1st order OSIP data.......' ttm=vis_osip(t,e[lm1].ccd_id,e[lm1].chipx,e[lm1].chipy,e[lm1].tg_mlam, order, fudge) ; colors ecol = 254 ; default ecol1 = 77 ; events, odd orders ecol2 = 122 ocol1 = 144 ; osip limits, odd order ocol2 = 188 ; osip limits, even order fcol1 = 199 ; osip fracresp color ecol21 = [ecol2, ecol1] ocol21 = [ocol2, ocol1] ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; IF (!ioflags.do_plot EQ 1) OR (!ioflags.do_plot EQ 99) THEN BEGIN ; plot 1st order MEG ; plot, e[lm].tg_mlam, $ abs(e[lm].tg_mlam) * e[lm].energy/!hc/1000., $ yra=[0,max(ords)+1],xtit='tg_mlam',ytit='order', title=fosip, col=ecol oplot, e[lm1].tg_mlam,$ abs(e[lm1].tg_mlam)*e[lm1].energy/!hc/1000.,$ col=ecol1 ; overplot the 1st order osip limits and fractional power ; oplot,e[lm1].tg_mlam,ttm.m_lo,col=ocol1 oplot,e[lm1].tg_mlam,ttm.m_hi,col=ocol1 oplot,e[lm1].tg_mlam,ttm.fracresp+(order-1),col=fcol1 ; loop over higher orders ; ll = where(ords NE 1, nn) ; already did 1st order FOR io = 0, nn-1 DO BEGIN order = ords[ll[io]] lmo=where( (abs(e.tg_m) EQ order) AND (e.tg_part EQ 2) AND lgood, n_lmo) ; IF n_lmo GE 1 THEN oplot, e[lmo].tg_mlam,abs(e[lmo].tg_mlam)*e[lmo].energy/!hc/1000.,col=88, psy = 4 message, /inf, 'Computing OSIP data for order='+string(order) ; use lm1, 1st order for boundary to ; get enough points; limits calc from ; m*lam for order m: ttm=vis_osip(t,e[lm1].ccd_id,e[lm1].chipx,e[lm1].chipy,e[lm1].tg_mlam, order, fudge) IF n_lmo GE 1 THEN plots,e[lmo].tg_mlam,abs(e[lmo].tg_mlam)*e[lmo].energy/!hc/1000.,col=ecol21[order MOD 2] oplot,e[lm1].tg_mlam,ttm.m_lo,col=ocol21[order MOD 2] oplot,e[lm1].tg_mlam,ttm.m_hi,col=ocol21[order MOD 2] oplot,e[lm1].tg_mlam,ttm.fracresp+(order-1),col=fcol1 ENDFOR ; save the plot: ; id_plot, 'MEG' IF !version.release GT 5.3 THEN BEGIN bmp_out, froot+'_MEG.bmp' ENDIF ELSE BEGIN gif_out, froot+'_MEG.gif' ENDELSE ENDIF IF (!ioflags.do_plot EQ 2) OR (!ioflags.do_plot EQ 99) THEN BEGIN ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Do same for HEG, starting at determining osip limits: ; message, /inf, 'Computing HEG 1st order OSIP data.......' order = 1 tth=vis_osip(t,e[lh1].ccd_id,e[lh1].chipx,e[lh1].chipy,e[lh1].tg_mlam, order, fudge) plot,e[lh].tg_mlam,abs(e[lh].tg_mlam)*e[lh].energy/!hc/1000.,yra=[0,max(ords)+1],xtit='tg_mlam',ytit='order', title=fosip, col = ecol oplot,e[lh1].tg_mlam,abs(e[lh1].tg_mlam)*e[lh1].energy/!hc/1000.,col=ecol1 ; overplot the osip limits and fractional probability: ; oplot,e[lh1].tg_mlam,tth.m_lo,col=ocol1 oplot,e[lh1].tg_mlam,tth.m_hi,col=ocol1 oplot,e[lh1].tg_mlam,tth.fracresp+(order-1),col=fcol1 ; loop over higher orders ll = where(ords NE 1, nn) ; already did 1st order FOR io = 0, nn-1 DO BEGIN order = ords[ll[io]] lho=where( (abs(e.tg_m) EQ order) AND (e.tg_part EQ 1) AND lgood, n_lho) ; IF n_lho GE 1 THEN plots, e[lho].tg_mlam,abs(e[lho].tg_mlam)*e[lho].energy/!hc/1000.,col=88,psy=4 message, /inf, 'Computing OSIP data for order='+string(order) ; use lh1, 1st order for boundary to ; get enough points; limits calc from ; m*lam for order m: tth=vis_osip(t,e[lh1].ccd_id,e[lh1].chipx,e[lh1].chipy,e[lh1].tg_mlam, order, fudge) IF n_lho GE 1 THEN oplot,e[lho].tg_mlam,abs(e[lho].tg_mlam)*e[lho].energy/!hc/1000.,col=ecol21[order MOD 2] oplot,e[lh1].tg_mlam,tth.m_lo,col=ocol21[order MOD 2] oplot,e[lh1].tg_mlam,tth.m_hi,col=ocol21[order MOD 2] oplot,e[lh1].tg_mlam,tth.fracresp+(order-1),col=fcol1 ENDFOR ; save the plot: ; id_plot, 'HEG' IF !version.release GT 5.3 THEN BEGIN bmp_out, froot+'HEG.bmp' ENDIF ELSE BEGIN gif_out, froot+'_HEG.gif' ENDELSE ENDIF ; close the ps file, if any: ; ps_close END ; E.g., to test: ; ;plt_osip,'acishf_001_osip-090_evt1a.fits','./Ard/acisD1999-07-22osipN0001.fits',lbl='Test',nevt=1000,orders=[1,2,3]