;% Time-stamp: <97/06/25 16:09:57 dph> ;% MIT Directory: ~dph/d3/XRCF_data/L1.5_devel ;% CfA Directory: /dev/null ;% File: rd.pro ;% Author: D. Huenemoerder ;% Original version: 970611 ;%===================================================================== ; NOTES: ; ; representing tg_part as bitmap w/ 3 bits per source ; limits to 30/3=10 sources. Could define wider bitmaps, ; but cumbersome. Use integer values? e.g., if could be ; HEG or MEG from same source, then instead of 0010 OR ; 0100, ??? ; I'll go to longs for src_map and part_map, limiting myself ; here to 10 sources. added src and part to refer to r,d ; assigned in the list. ; ;;; actually - don't need src_map if src encoded into part... ; Details of mapping is implementation issue, rather than ; logical. ; 970613 - long talks w/ Dan D. about L1.5 file, collisions - remove ; ambiguity if possible, otherwise, leave it and don't provide r,d - ; NO overlap extensions. ; Also - make more compatible w/ LETG/HRC via m=+1, -1, when ; un-resolved; could use as HEG+MEG ; ; IDL main to do L1.5 prototype. ;;;;;;;;;;;;;;;;;;;;;;;;;; IF n_elements(do_read) EQ 0 THEN do_read = 1 IF n_elements(do_plot) EQ 0 THEN do_plot = 99 IF n_elements(do_ps) EQ 0 THEN do_ps = 0 IF n_elements(do_gif) EQ 0 THEN do_gif = 0 ;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;; ; Use a smallish molecular ; contamination ACIS-S L1 file: ;fname = 'acis_hhasmc20001i002n000_evt.fits' fname = 'acis-s_evt.fits' ev_ext = 2 ; the event extension ; read the events into a structure: IF do_read THEN BEGIN ; data-preperation ev = mrdfits(fname, ev_ext, ev_hdr) ev_tags = tag_names(ev) nx_tag = (where(ev_tags EQ 'X')+1)[0] ; force to scalar ny_tag = (where(ev_tags EQ 'Y')+1)[0] x_scale = sxpar(ev_hdr, 'TCDLT'+ strcompress(string(nx_tag), /rem)) y_scale = sxpar(ev_hdr, 'TCDLT'+ strcompress(string(ny_tag), /rem)) ;IDL> help,/st,ev ;** Structure <37fd90>, 15 tags, length=48, refs=1: ; TIME DOUBLE 1.0463936e+08 ; EXPNO LONG 316 ; CCD_ID INT 6 ; CHIPX INT 51 ; CHIPY INT 1022 ; TDETX INT 11473 ; TDETY INT 13618 ; DETX INT 5265 ; DETY INT 3797 ; X INT 5266 ; Y INT 3797 ; PHA LONG 235 ; PI LONG 904 ; GRADE INT 0 ; FLTGRADE INT 0 nil = mrdfits(fname, 0, prim_hdr) ; get the primary header ;;;;;;;;;;;;;;;;;;;;;;;;;; ; truncate the data to smaller size, ; only for prototyping ease n_el = 30000 ; a reasonable number, based on inspection. i_time = sort(ev.time) ; time-sort: ev.time = ev[i_time].time ev.EXPNO = ev[i_time].EXPNO ev.CCD_ID = ev[i_time].CCD_ID ev.CHIPX = ev[i_time].CHIPX ev.CHIPY = ev[i_time].CHIPY ev.TDETX = ev[i_time].TDETX ev.TDETY = ev[i_time].TDETY ev.DETX = ev[i_time].DETX ev.DETY = ev[i_time].DETY ev.X = ev[i_time].X ev.Y = ev[i_time].Y ev.PHA = ev[i_time].PHA ev.PI = ev[i_time].PI ev.GRADE = ev[i_time].GRADE ev.FLTGRADE = ev[i_time].fltgrade ev = ev[0:n_el-1] ; truncate ;;;;;;;;;;;;;;;;;;;;;;;;;; ; STUB ; Source "DETECTION" happens ; here..... get the source ; coordinates. Since these aren't ; quite right in the xrcf header, fix ; them by adding a priori known ; offset. foax,foay is zero-order ; coordinate in the "sky" coordinate ; system. ; ; For flight - source is not ; necessarily at foax, foay. ; Detection will produce a list of ; sources and coordinates in sky x,y. foax = sxpar(prim_hdr, 'FOAFP_X') foay = sxpar(prim_hdr, 'FOAFP_Y') foay = 4237.5 ; hack; empirically determined. x_src = foax y_src = foay ;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;; ; STUB: read aspect solution, ; The aspect solution will give ; a series of vectors: time, ; delta_ra, delta_dec, ; delta_roll, and uncertainties ; on each. The delta- is ; relative to requested (I ; think; but definitely not ; (ra,dec)=(0,0), since that is ; a singular system at the pole) ; ; determine mean roll angle ; relative to axaf_y, and range ; of roll. (Is this ROLL_NOM ; from the header? Probably not ; - that is the intended value, ; not determined from aspect.) ; mean_roll_angle = avg(roll) ; in practise, more like this ; mean_roll_angle = 2. ; degrees ccw from det_x == axaf_y mean_roll_angle = 5. ; degrees ccw from det_x == axaf_y ; (set it to something small here for ; testing) ; min_roll=min(roll) ; max_roll=max(roll) ; roll_range=max_roll-min_roll ; need to check range to see if out of ; limits for region definition. roll_range = 0.5 ; degrees about mean_roll_angle (half-cone) ; recent memo from Tom Aldcroft ; indicates that roll-drift in 48 ; hours will be only about 0.05 ; degrees (3 arcmin). ;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;; ; STUB: alter the coordinates a bit, ; so they are more like sky ; coords. assume roll is about ; source, which isn't ; necessarily the case. ; change the sky coords a bit - rotate ; x,y. dx = randomu(seed, n_el) ; anti-aliasing terms; randomize ; within a pixel. dy = randomu(seed, n_el) rotxy, ev.x-x_src+dx, ev.y-y_src+dy, mean_roll_angle, x, y x = x+x_src y = y+y_src ;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;; ; turn x,y into angles from aim-point x = (x-x_src)*double(x_scale) * !dtor ; radians, from src y = (y-y_src)*double(y_scale) * !dtor ; radians ; the event data is now more ; flight-like - no dither, but x,y not ; aligned w/ detx, dety, chipx, chipy. ; ALSO: assuming x,y have been changed ; to reals and angles. This is ; TBD. May at least be longs. ; define an output structure: evt = {event, $ time: 1.d0,$ EXPNO: 1L, $ CCD_ID: 1, $ CHIPX: 1, $ CHIPY: 1, $ TDETX: 1, $ TDETY: 1, $ DETX: 1, $ DETY: 1, $ X: 1.0, $ Y: 1.0, $ PHA: 1L, $ PI: 1L, $ GRADE: 1, $ FLTGRADE:1, $ ; tg_src_map: 0L, $ ; L1.5 additions start here --------- tg_part_map: 0L, $ tg_r: 0.0, $ tg_d: 0.0, $ tg_mlam: 0.0, $ tg_order: 0, $ tg_lam:0.0 $ } event_list = replicate( {event}, n_el) ; copy input to output event_list.time = ev.time ; copy event_list.expno = ev.expno ; copy event_list.ccd_id = ev.ccd_id ; copy event_list.chipx = ev.chipx ; copy event_list.chipy = ev.chipy ; copy event_list.tdetx = ev.tdetx ; copy event_list.tdety = ev.tdety ; copy event_list.detx = ev.detx ; copy event_list.dety = ev.dety ; copy event_list.x = x ; new float angle event_list.y = y ; new float angle event_list.pha = ev.pha ; copy event_list.pi = ev.pi ; copy event_list.grade = ev.grade ; copy event_list.fltgrade = ev.fltgrade ; copy ; ; here on get filled in at Level 1.5 ; event_list.tg_src_map = 0 ; no sources event_list.tg_part_map = 0 ; no parts yet, all bkg. event_list.tg_r = 999. ; bogus value event_list.tg_d = 999. ; bogus value event_list.tg_mlam = 999. ; bogus value event_list.tg_order = 255 ; bogus value event_list.tg_lam = 0.0 ; bogus value ; we can now remove the input ; structure and work on the output: ; punt unused vectors to conserve memory ev = 0B x = 0B y = 0B dx = 0B dy = 0B i_time = 0B do_read = 0 ; data has been read and munged. ENDIF ; END DATA-PREPARATION ;;;;;;;;;;;;;;;;;;;;;;;;;; ; a clean slate: tag_non_diffracted = 999 ; tag for no order assignment (999 is ; used by HETG) ;event_list.tg_src_map = 0 ; no sources event_list.tg_part_map = 0 ; no parts yet, all bkg. event_list.tg_r = 999. ; bogus value event_list.tg_d = 999. ; bogus value event_list.tg_mlam = 999. ; bogus value event_list.tg_order = tag_non_diffracted ; bogus value (0, -1 are valid values!) event_list.tg_lam = 0.0 ; bogus value ;;;;;;;;;;;;;;;;;;;;;;;;;; ; L1.5 processing begins here. ;;;;;;;;;;;;;;;;;;;;;;;;;; ; make up src bit fields within tg_part src_field = '0007'XL src_factor = 8L src_1_bits = src_field src_2_bits = src_1_bits * src_factor ; bit-shift left by field-width src_3_bits = src_2_bits * src_factor ; ...and again for next src. ; etc ; STUB: determine regions sizes for ; HEG, MEG. Depends on several ; factors. ; astigmatic width: R_TG_MEG = 600. ; approx mm radius MEG outer ring ; (this will come from reference data) R_TG_HEG = 425. ; approx mm radius HEG outer ring x_max = max(abs(event_list.x)) ; approximate maximum for tg_r foc_len = sxpar(prim_hdr,'FOC_LEN') ; imaging focal length, [mm]; this ; is in header, but could also ; come from ref data. tg_foc_len = 8782.8d0 ; should be in header, but not.; also ; come from ref data. dfocus = sxpar(prim_hdr, 'DEFOCUS') ; mm defocus, for particular ; observation, should be in header. ; calculate cross-dispersion width due ; to astigmatism and de-focus: term = 2.d0 / tg_foc_len * ( x_max^2 + abs(dfocus)/tg_foc_len ) d_astig_dfoc_meg = term * R_TG_MEG ; radians d_astig_dfoc_heg = term * R_TG_HEG ; radians ; calculate cross-dispersion width (in ; sky coords) due to roll range ; (probably very small, in practice, ; but must allow for full range in ; specs). d_roll = roll_range * !dtor * x_max ; STUB: ; source extent; approx point source ; psf as a couple arcsec; could come ; from source detection ouput. hrma_fwhm = 1.d0 ; arcsec, approx fwhm d_hrma_psf = hrma_fwhm /3600. * !dtor ; radians ; combine results box_factor = 2.*sqrt(3) ; fwhm-to-sigma factors gauss_factor = 2.*sqrt(2.*alog(2)) width_factor = 6. ; how many sigmas to use. sig_cross_disp_meg = width_factor * $ sqrt( $ (d_hrma_psf/gauss_factor)^2 + $ (d_astig_dfoc_meg/box_factor)^2 $ ) sig_cross_disp_heg = width_factor * $ sqrt( $ (d_hrma_psf/gauss_factor)^2 + $ (d_astig_dfoc_heg/box_factor)^2 $ ) meg_d_width = sig_cross_disp_meg + d_roll heg_d_width = sig_cross_disp_heg + d_roll ;;;;;;;;;;;;;;;;;;;;;;;;;; ; STUB ; get grating angles !meg_angle = 4.65d0 ; degrees clocking angle. this is ; both grating and detector combined. ; should be kept separate; will come ; from ref data for grating and ; detector, be in pixlib params. !heg_angle = -5.35d0 ; rotate x,y to meg and heg approx r,d ; planes, for selection of photons via ; sky x,y. rotxy, event_list.x, event_list.y, $ !meg_angle-mean_roll_angle, r_meg, d_meg rotxy, event_list.x, event_list.y, $ !heg_angle-mean_roll_angle, r_heg, d_heg r_zero_ord = sqrt(event_list.x^2+event_list.y^2) ; photon angle from zero-order ; centroid. ;;;;;;;;;;;;;;;;;;;;;;;;;; ; assign region bitmap values ; select regions ;event_list.tg_src_map = 0 ; Only one source in this test. ; If another, this is a bitmap showing ; which source a photon could be ; from. tag_bkg = 0B tag_zero_order = 1B tag_meg = 2B tag_heg = 4B ; additional sources in field would ; have 8*(srcnum-1) added to these ; values. (if sources started ; enumeration at 1) This ; implementation uses type LONG for ; the bitmap, so limits us to max of 8 ; sources. l_meg = where( abs(d_meg) LT meg_d_width) l_heg = where( abs(d_heg) LT heg_d_width) n_hrma_psf = 50 ; should be big enough! (but make param) l_zero = where( r_zero_ord LT d_hrma_psf * n_hrma_psf) r_zero_ord = 0B ; punt; recover memory ; apply region logic ; meg and heg are allowed to overlap. ; Zero-order has precedence. event_list[l_meg].tg_part_map = $ event_list[l_meg].tg_part_map OR tag_meg event_list[l_heg].tg_part_map = $ event_list[l_heg].tg_part_map OR tag_heg event_list[l_zero].tg_part_map = tag_zero_order ;;;;;;;;;;;;;;;;;;;;;;;;;; ; Find overlap regions: ; For one source only, collision ; between HEG and MEG is: l_overlap = where( event_list.tg_part_map EQ (tag_meg OR tag_heg)) l_meg = where(event_list.tg_part_map EQ tag_meg) ; re-define meg only part l_heg = where(event_list.tg_part_map EQ tag_heg) ; re-define heg only part ;;;;;;;;;;;;;;;;;;;;;;;;;; IF (do_plot EQ 1) OR (do_plot EQ 99) THEN begin ; let's see it! wshow col_meg = 135 col_heg = 100 col_lap = 190 col_bkg = !d.n_colors-1 col_zero = 80 plot, event_list.x*!radeg, event_list.y*!radeg, $ xra = [-2,2]/1000.*!radeg, yra = [-2,2]/2000.*!radeg, /xst, /yst, $ xtit = 'Sky X [degrees]', ytit = 'Sky Y [degrees', $ title = 'Geometric identification of spectral parts' oplot, event_list[l_meg].x*!radeg, event_list[l_meg].y*!radeg, col = col_meg oplot, event_list[l_heg].x*!radeg, event_list[l_heg].y*!radeg, col = col_heg oplot, event_list[l_zero].x*!radeg, event_list[l_zero].y*!radeg, col=col_zero oplot, event_list[l_overlap].x*!radeg, event_list[l_overlap].y*!radeg, $ psy = 1, col = col_lap xyouts, 0.1, 0.005, 'MEG+', col = col_meg, ali = 0.5 xyouts, 0.1, -0.03, 'HEG+', col = col_heg, ali = 0.5 xyouts, -0.1, 0.02, 'HEG-', col = col_heg xyouts, -0.1, -0.005, 'MEG-', col = col_meg xyouts, 0, 0.02, 'Zero order', align = 0.5, col = 80 xyouts, 0.06, -0.0075, 'overlapped', col = col_lap, chars = 0.8 xyouts, 0.01, -0.04, 'background' id_plot, 'dph' IF do_gif THEN wgif, !window, 'Sky_xy_regions.gif' ENDIF ;;;;;;;;;;;;;;;;;;;;;;;;;; ; determine scale factor: ; Need to convert detxy to angular ; units (use radians for now) Scale ; isn't in header - should be in a ; TCDLT keyword of the events ; extension. ndx_tag = (where(ev_tags EQ 'DETX')+1)[0] ; force to scalar ndy_tag = (where(ev_tags EQ 'DETY')+1)[0] dx_scale = sxpar(ev_hdr, 'TCDLT'+ strcompress(string(ndx_tag), /rem)) dy_scale = sxpar(ev_hdr, 'TCDLT'+ strcompress(string(ndy_tag), /rem)) ; if dx_scale isn't the same as ; dy_scale, then we have a serious ; problem (w/ the FITS header) dx_scale = dx_scale * !dtor ; degrees/pixel to radians/pixel pixsize = 0.024 ; mm dx_scale = pixsize/tg_foc_len ; STUB; radians/pixel ;;;;;;;;;;;;;;;;;;;;;;;;;; ; Now compute r,d. Photons tagged by ; part of spectrum. Now use chipid, ; chipx, chipy to compute r,d w/ ; pixlib calls. ; STUB - ; do in detx,dety by rotating by ; empirical angles (instead of ; letting pixlib to detailed ; tranformations using full ; geometry). rotxy, event_list[l_meg].detx-x_src, $ event_list[l_meg].dety-y_src, $ !meg_angle, r_meg, d_meg rotxy, event_list[l_heg].detx-x_src, $ event_list[l_heg].dety-y_src, $ !heg_angle, r_heg, d_heg ; combine, r, d into one vector event_list[l_meg].tg_r = r_meg * dx_scale event_list[l_meg].tg_d = d_meg * dx_scale event_list[l_heg].tg_r = r_heg * dx_scale event_list[l_heg].tg_d = d_heg * dx_scale ; assign overlap coords to MEG, since ; MEG spectrum comes in closer. This ; is arbitrary choice. Order-sorting ; per arm will re-identify. rotxy, event_list[l_overlap].detx-x_src, $ event_list[l_overlap].dety-y_src, $ !meg_angle, r_meg, d_meg event_list[l_overlap].tg_r = r_meg * dx_scale event_list[l_overlap].tg_d = d_meg * dx_scale ; Compute alternate coords for ; overlapped events. event_list_overlap = event_list[l_overlap] ; copy to new structure rotxy, event_list_overlap.detx-x_src, $ event_list_overlap.dety-y_src, $ !heg_angle, r_heg, d_heg event_list_overlap.tg_r = r_heg * dx_scale event_list_overlap.tg_d = d_heg * dx_scale ;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;; IF (do_plot EQ 2) OR (do_plot EQ 99) THEN begin ; show it all wset,1 wshow sav_pmulti = !p.multi !p.multi = [1, 1, 2] erase ; im = make_image(event_list.x*!radeg, event_list.y*!radeg, $ ; xra = [-0.5, .5], yra = [-.5, .5]/2,$ ; xbin = 0.001, ybin = 0.001, $ ; xax = axx, yax = axy) plot, event_list.x*!radeg, event_list.y*!radeg, $ xra = [-0.5, .5], $ ; yra = [-.1, .1], /xst, /yst, $ xtit = 'Sky X [degrees]', ytit = 'Sky Y [degrees]', $ title = 'Field View: HETG/ACIS-S' ll = where( event_list.tg_part_map EQ tag_meg ) oplot, event_list[ll].x*!radeg, event_list[ll].y*!radeg;, col = 180 ll = where( event_list.tg_part_map EQ tag_heg ) oplot, event_list[ll].x*!radeg, event_list[ll].y*!radeg;, col = 100 ll = where( event_list.tg_part_map AND tag_zero_order ) oplot, event_list[ll].x*!radeg, event_list[ll].y*!radeg;, col = 80 ; ll = where( (event_list.tg_part_map EQ tag_bkg) ) ; oplot, event_list[ll].x, event_list[ll].y, col = 55 ll = where( event_list.tg_part_map EQ (tag_meg OR tag_heg) ) oplot, event_list[ll].x*!radeg, event_list[ll].y*!radeg, psy = 1;col = 144 id_plot, 'dph' IF do_gif THEN wgif, !window, 'Field_view.gif' wset,0 !p.multi = sav_pmulti ENDIF ;;;;;;;;;;;;;;;;;;;;;;;;;; ; STUB: ; compute m*lambda. Note precedence ; order, for overlapped regions isn't ; treated explicitly. The second ; assignment of tg_mlam will overwrite ; some of the first, if there are ; overlapped regions. ; ; Use simple grating equation: ; m*lambda = Period * sin(beta) ; ~ Period * tg_r[mm] / X_TG ; = Period * tg_r ; Period = meg or heg period ; tg_r[mm] = diffraction coord in mm ; tg_r = diffraction coord in radians ; X_TG = grating focal length ; (Rowland diameter) ; ; (This will really be done via pixlib ; calls which treat geometry in detail.) event_list.tg_mlam = event_list.tg_r event_list[l_heg].tg_mlam = !heg_period * event_list[l_heg].tg_r event_list[l_meg].tg_mlam = !meg_period * event_list[l_meg].tg_r ; ditto for overlap region, except a ; priori determined to be all HEG: event_list_overlap.tg_mlam = !heg_period * event_list_overlap.tg_r ;;;;;;;;;;;;;;;;;;;;;;;;;; ; Now it's time for order-sorting, if ; ACIS. ; An order estimator is: ; m_est = PI/hc * (m*lambda). ; ; A tolerance about the estimate must ; be specified. This will come from a ; PI RMF (response matrix, in PI ; units). Here I will use a ballpark ; CCD resolution: E/dE ~ 50 ~ m/dm. ; So, dm = m/50 (one-sigma). To make ; wider, multiply by some factor. ; STUB ; ccd_fwhm is really ccd_fwhm(E). ; Could assume sqrt law here ; (theoretical value) - or MARX FI/BI ; functions. I practice, will come ; from an RMF or other calibration data. ; next two lines are for constant E/dE ccd_fwhm = 50. ; E/dE = m/dm, fwhm ccd_sigma = ccd_fwhm/gauss_factor ; convert to sigma ; If a function of energy is used, ; then the calculation is a bit more ; complicated because the ccd_sigma to ; apply then depends on the order. n_sigma_hi = 4.5 ; How wide are regions in PI? These ; are also energy-dependent. They ; will be determined from the RMF (or ; equivalent data) in terms of x% ; enclosed area. Since the ; distribution is asymmetric, the high ; and low limits can be different ; (when centered on the maximum or ; gaussian-fit centroid). Care will ; have to be taken to deal with ; "collisions": a low-limit ; overlapping with the high-limit of ; the order below. n_sigma_lo = 4.5 m_estimator = (event_list.pi/1000.) / !hc * event_list.tg_mlam ; pi is in eV. max_order = 10 ; max absolute value of order to ; process; should be parameter. event_list.tg_order = tag_non_diffracted ; no order assigned yet event_list[l_zero].tg_order = 0 ; assign zero-order tag ; ditto for overlap region ; xxx FLAWED LOGIC HERE: ; m_guess_overlap is wrong - it could ; be either HEG or MEG, so need a pair ; of order estimators! and propagate ; affect onward... m_estimator_overlap = (event_list_overlap.pi/1000.) / !hc * $ event_list_overlap.tg_mlam event_list_overlap.tg_order = tag_non_diffracted ; For energy dependent ccd_sigma, use ; m_estimator and m*lambda to estimate ; photon E and look up sigma(E). ; ; STUB: use simple model, instead of ; real ccd_sigma(E) energy = findgen(1000)/1000.*12000. ; energy vector, 0 to 12000 eV ; Simple CCD noise model: ccd_gain = 3.68 ; eV / electron ccd_fano = 0.115 ; fano factor ccd_rnoise = 4.0 ; read_noise, electrons ccd_sigma_model = ccd_gain * sqrt( ccd_rnoise^2 + ccd_fano*energy/ccd_gain) ccd_sigma_model = energy / ccd_sigma_model ; convert to E/dE. ; now "smooth" m_estimator and use to ; calc a ccd_sigma per photon: ; interpolate approx energy from model: m_guess = round(m_estimator) m_guess_overlap = round(m_estimator_overlap) ccd_sigma = interpol_sort(ccd_sigma_model, energy, $ (1000.*!hc/event_list.tg_mlam) * m_guess) ccd_sigma_overlap = interpol_sort(ccd_sigma_model, energy, $ (1000.*!hc/event_list_overlap.tg_mlam) * $ m_guess_overlap ) ; compute upper and lower limits to ; ratios of m_estimator/m_integer upper_ratio = 1. + n_sigma_hi / ccd_sigma lower_ratio = 1. - n_sigma_lo / ccd_sigma ; calculate actual ratios obs_ratio = m_estimator / m_guess ; If the obs_ratio is between the ; lower_ratio and upper_ratio, then ; the order is identified. (BUT - ; there is no check for collision). l_identified = where( (obs_ratio LT upper_ratio) AND $ (obs_ratio GE lower_ratio) ) IF (size(l_identified))[0] THEN $ ; there are some event_list[l_identified].tg_order = m_guess[l_identified] ; fix 0/0 points - replace w/ zero again: event_list[l_zero].tg_order = 0 ;;;;;;;;;;;;;;;;;;;; ; ditto for overlap... ; compute upper and lower limits to ; ratios of m_estimator/m_integer ; save the zero-order indices: l_zero_ov = where(event_list_overlap.tg_part_map) EQ tag_zero_order upper_ratio_ov = 1. + n_sigma_hi / ccd_sigma_overlap lower_ratio_ov = 1. - n_sigma_lo / ccd_sigma_overlap ; calculate actual ratios obs_ratio_ov = m_estimator_overlap / m_guess_overlap ; If the obs_ratio is between the ; lower_ratio and upper_ratio, then ; the order is identified. l_identified_ov = where( (obs_ratio_ov LT upper_ratio_ov) AND $ (obs_ratio_ov GE lower_ratio_ov) ) IF (size(l_identified_ov))[0] THEN $ ; there are some event_list_overlap[l_identified_ov].tg_order = $ m_guess_overlap[l_identified_ov] ; replace zeroes... IF (size(l_zero_ov))[0] THEN $ ; there are some event_list_overlap[l_zero_ov].tg_order = 0 ;;;;;;;;;;;;;;;;;;;;;;;;;; ; see it IF (do_plot EQ 3) OR (do_plot EQ 99) THEN BEGIN wshow plot, event_list.tg_mlam, m_estimator, yra = [-10,10], xra = [-35,35], /nod plots, event_list.tg_mlam, m_estimator, color = $ abs(event_list.tg_order)*20, noclip = 0 ENDIF ;;;;;;;;;;;;;;;;;;;;;;;;;; ; and at last, wavelengths! ; look out for /0, 0/0, or other ; undefined points: skip zero-order, ; background in ratio. lmh = where( (event_list.tg_order NE tag_non_diffracted) $ AND $ (event_list.tg_order NE 0) $ ) event_list[lmh].tg_lam = event_list[lmh].tg_mlam / event_list[lmh].tg_order lmh = where( (event_list_overlap.tg_order NE tag_non_diffracted) $ AND $ (event_list_overlap.tg_order NE 0) $ ) event_list_overlap[lmh].tg_lam = $ event_list_overlap[ lmh ].tg_mlam / event_list_overlap[ lmh ].tg_order ;;;;;;;;; memory cleanup... m_guess = 0B m_guess_overlap = 0B upper_ratio = 0B lower_ratio = 0B obs_ratio = 0B l_identified = 0B upper_ratio_ov = 0B lower_ratio_ov = 0B obs_ratio_ov = 0B l_identified_ov = 0B m_estimator = 0B l_meg = 0B l_heg = 0B energy = 0B ccd_sigma = 0B ;;;;;;;;;;;;;;;;;;;;;;;;;; ; Need to re-assign overlap region ; photons. Those with an order ; assigned only as MEG or only as HEG ; may be considered unambiguous ; identifications and will have all ; coordinates specified. Those which ; are ambiguous will have multiple field ; bits set in tg_part_map and ; which indicate the src, photon ; combinations ; Indeterminate points will not have ; r, d, m*lambda, m, or lambda set to ; valid values. (m=999, the rest 0.0 ; TBR) ; ; l_overlap is the array of ; indices of overlap region photons in ; the original list. ; ; Compare: ; ; a = event_list[l_overlap].tg_order ; ; b = event_list_overlap.tg_order ; ; Case 1) ; IF a NE tag_non_diffracted ; AND ; b EQ tag_non_diffracted then ; identified w/ a: clear b's bits in ; tg_part_map. ; ; Case 2) ; IF a EQ tag_non_diffracted ; AND ; b EQ tag_non_diffracted ; background to both. ; Clear tg_part_map ; ; Case 3) ; IF a EQ tag_non_diffracted ; AND ; b NE tag_non_diffracted, then ; identified w/ b; ; swap a,b (all params, not just ; t[g_order - though non-L1.5 ; columns should be the same if ; indexed correctly!); ; clear a's bits in tg_part_map. ; ; Case 4) ; IF a NE tag_non_diffracted ; AND ; b NE tag_non_diffracted, then ; ambiguous; leave tg_part_map, ; tg_src_map pointing to both; clear ; r,d, etc coords. ; ; Bit clearing/setting- careful to ; handle same-src parts (e.g., MEG ; collides w/ HEG from same source). ; This example only treats one ; source, so deal only w/ ; tg_part_map first field (bits 111). message, 'Process overlap: # events '+string(n_elements(l_overlap)), /inf a_tmp = event_list[l_overlap].tg_order b_tmp = event_list_overlap.tg_order ; case 1: all from MEG part: ll = where( (a_tmp NE tag_non_diffracted) $ AND $ (b_tmp EQ tag_non_diffracted) ) message, 'Overlap Case 1: # only MEG '+string(n_elements(ll)), /inf l_1 = l_overlap[ll] ; save for plotting ; NOTE- tricky indexing. Have ; truncated event_list to l_overlap ; elements, then truncate again. Have ; to double-index to get them back in ; the right place. IF (size(ll))[0] THEN $ ; there are some... event_list[l_1].tg_part_map = tag_meg ; overlap was in HEG coords ; and None were identified as HEG. ; case 2: both background ll = where( (a_tmp EQ tag_non_diffracted) $ AND (b_tmp EQ tag_non_diffracted) ) message, 'Overlap Case 2: # both bkg '+string(n_elements(ll)), /inf l_2 = l_overlap[ll] IF (size(ll))[0] THEN $ ; there are some... event_list[l_2].tg_part_map = tag_bkg ; case 3: Identified as HEG - swap: ll = where( (a_tmp EQ tag_non_diffracted) $ AND (b_tmp NE tag_non_diffracted) ) message, 'Overlap Case 3: # only HEG '+string(n_elements(ll)), /inf l_3 = l_overlap[ll] IF (size(ll))[0] THEN BEGIN ; there are some... event_list[l_3] = event_list_overlap[ll] event_list[l_3].tg_part_map = tag_heg ENDIF ; Case 4: ambiguous ll = where( (a_tmp NE tag_non_diffracted) $ AND (b_tmp NE tag_non_diffracted) ) message, 'Overlap Case 4: # ambiguous '+string(n_elements(ll)), /inf l_4 = l_overlap[ll] IF (size(ll))[0] THEN BEGIN ; there are some... event_list[l_4].tg_r = 999. event_list[l_4].tg_d = 999. event_list[l_4].tg_mlam = 999. event_list[l_4].tg_order = tag_non_diffracted event_list[l_4].tg_lam = 999. ENDIF message, 'l_1: '+ string(n_elements(l_1)) + ' + ', /inf ; debug message, 'l_2: '+ string(n_elements(l_2)) + ' * ', /inf ; debug message, 'l_3: '+ string(n_elements(l_3)) + ' dia ', /inf ; debug message, 'l_4: '+ string(n_elements(l_4)) + ' tri ', /inf ; debug ;;;;;;;;;;;;;;;;;;;;;; IF (do_plot EQ 4) OR (do_plot EQ 99) THEN begin ; visual check surface, dist(20), xra=[-1,1.2]*800, $ yra=[-1,1]*300, $ zra=[0,8], /nod, /sav, chars = 1.6, ax = 70, az = -20 llh = where( (abs(event_list.tg_order) GT 0) AND $ (abs(event_list.tg_order) LE 10) AND $ event_list.tg_part_map EQ tag_heg) llm = where( (abs(event_list.tg_order) GT 0) AND $ (abs(event_list.tg_order) LE 10) AND $ event_list.tg_part_map EQ tag_meg) xm = event_list[llm].detx-x_src ym = event_list[llm].dety-y_src xh = event_list[llh].detx-x_src yh = event_list[llh].dety-y_src plots, xm, ym, noclip=0, /t3d, col=66 plots, xh, yh, noclip=0, /t3d, col=166 ch = 166 ;abs(event_list[llh].tg_order) * 80 cm = 66 ;abs(event_list[llm].tg_order) * 20 plots, xm, ym, event_list[llm].pi/1000., noclip=0, /t3d, col=cm plots, xh, yh, event_list[llh].pi/1000., noclip=0, /t3d, col=ch plots, event_list[l_1].detx-x_src, $ event_list[l_1].dety-y_src, $ event_list[l_1].pi/1000., $ noclip = 0, /t3d, col = 144, psy = 1, symsi = 2 plots, event_list[l_2].detx-x_src, $ event_list[l_2].dety-y_src, $ event_list[l_2].pi/1000., $ noclip = 0, /t3d, col = 111, psy = 2, symsi = 2 plots, event_list[l_3].detx-x_src, $ event_list[l_3].dety-y_src, $ event_list[l_3].pi/1000., $ noclip = 0, /t3d, col = 88, psy = 4, symsi = 2 plots, event_list[l_4].detx-x_src, $ event_list[l_4].dety-y_src, $ event_list[l_4].pi/1000., $ noclip = 0, /t3d, col = 66, psy = 5, symsi = 2 ENDIF ;;;;;;;;;;;;;;;;;;;;;;;; ; make some more figures ; r,d plot: IF (do_plot EQ 5) OR (do_plot EQ 99) THEN BEGIN sav_pmulti = !p.multi !p.multi = [0, 1, 2] l_meg = where(event_list.tg_part_map EQ tag_meg) l_heg = where(event_list.tg_part_map EQ tag_heg) plot, event_list[l_meg].tg_r*!radeg, $ event_list[l_meg].tg_d*!radeg, $ xra = [-.5, .5], yra = [-.005, 0.005], $ ; xtit = 'tg_r: dispersion angle [degrees]', $ ; ytit = 'tg_d: cross-dispersion angle !C[degrees]', $ title = 'MEG events in diffraction coordinates' plot, event_list[l_heg].tg_r*!radeg, $ event_list[l_heg].tg_d*!radeg, $ xra = [-.5, .5], yra = [-.005, 0.005], $ xtit = 'tg_r: dispersion angle [degrees]', $ ytit = ' tg_d: cross-dispersion angle [degrees]', $ title = 'HEG events in diffraction coordinates' id_plot, 'dph' IF do_gif THEN wgif, !window, 'Diffraction_coords.gif' !p.multi = sav_pmulti ENDIF ; m*lambda plot IF (do_plot EQ 6) OR (do_plot EQ 99) THEN BEGIN sav_pmulti = !p.multi !p.multi = [0, 1, 2] l_meg = where(event_list.tg_part_map EQ tag_meg) l_heg = where(event_list.tg_part_map EQ tag_heg) plot, event_list[l_meg].tg_mlam, $ event_list[l_meg].tg_d*!radeg, $ xra = [-30, 30], yra = [-.005, 0.005], /xst, /yst, $ title = 'MEG events in first-order wavelengths' plot, event_list[l_heg].tg_mlam, $ event_list[l_heg].tg_d*!radeg, $ xra = [-15, 15], yra = [-.005, 0.005], /xst, /yst, $ xtit = 'tg_mlam: m*lambda ['+!Angstrom+']', $ ytit = ' tg_d: cross-dispersion angle [degrees]', $ title = 'HEG events in first-order wavelengths' id_plot, 'dph' IF do_gif THEN wgif, !window, 'Wavelength_m1_coords.gif' !p.multi = sav_pmulti ENDIF ; pi plot. IF (do_plot EQ 7) OR (do_plot EQ 99) THEN BEGIN ; sav_pmulti = !p.multi ; !p.multi = [0, 1, 2] l_meg = where( (event_list.tg_part_map EQ tag_meg) AND $ (event_list.tg_order NE tag_non_diffracted)) l_heg = where( (event_list.tg_part_map EQ tag_heg) AND $ (event_list.tg_order NE tag_non_diffracted)) l_nil = where( (event_list.tg_part_map EQ tag_meg) $ AND $ (event_list.tg_order EQ tag_non_diffracted)) plot, event_list[l_meg].tg_mlam, $ event_list[l_meg].pi/1000., $ xra = [0, 30], yra = [0, 10], /xst, /yst, $ xtit = 'tg_mlam: m*lambda ['+!angstrom+']', $ ytit = 'CCD PI [keV]', $ title = 'Order sorting MEG events...', /nod mcol = !d.n_colors - (event_list[l_meg].tg_order)*20 plots, event_list[l_meg].tg_mlam, $ event_list[l_meg].pi/1000., $ color = mcol, $ noclip = 0, psy = 3 plots, event_list[l_nil].tg_mlam, $ event_list[l_nil].pi/1000., psy = 3, noclip = 0 id_plot, 'dph' IF do_gif THEN wgif, !window, 'PI_mlam_meg.gif' ; !p.multi = sav_pmulti ENDIF ; counts plot IF (do_plot EQ 8) OR (do_plot EQ 99) THEN BEGIN lm1=where( (event_list.tg_part_map eq tag_meg) $ and $ (abs(event_list.tg_order) eq 1)) lh1=where( (event_list.tg_part_map eq tag_heg) $ and $ (abs(event_list.tg_order) eq 1)) sm1=reform( make_image(event_list[lm1].tg_lam, event_list[lm1].tg_d,$ xra=[2,25], yra=[-1,1]*2.e-5,$ xbin=0.05,ybin=4.e-5,xax=axm) ) sh1=reform( make_image(event_list[lh1].tg_lam, event_list[lh1].tg_d,$ xra=[2,25]/2, yra=[-1,1]*2.e-5,$ xbin=0.05,ybin=4.e-5,xax=axh) ) plot,axm,sm1,psy=10,/yty,yra=[1,100],xra=[0,25], /xst, $ xtit = 'Wavelength ['+!angstrom+']', $ ytit = 'Counts per bin [0.05 '+!angstrom+' bins]', $ tit = 'MEG and HEG first order counts' oplot,axh,sh1,psy=10,col=col_heg xyouts, 20, 20, 'MEG m=+-1';, col = col_meg xyouts, 5, 50, 'HEG m=+-1', col = col_heg id_plot, 'dph' IF do_gif THEN wgif, !window, 'Cts_spec.gif' ENDIF ;;;;;;;;;;;;;;;;;;;;;;;;;; ; finally, write a FITS bintable. ; Hack this - use other tools to fix ; header (e.g., ftools, fmodhead, ; etc), and copy GTI extension. ;;;;;;;;;;;;;;;;;;;;;;;; END