function pha_grp_histo, data ndata = n_elements(data) n = where(data.grouping eq 1, ngroup) nhisto = 2*ngroup histo = replicate({wave:0D, flux:0D, grp_num:0}, nhisto) for i=0, ngroup-1 do begin first_row = n[i] last_row = n[i]+data[n[i]].chans_per_grp-1 ;print, row_low, row_high histo[i*2].wave = data[first_row].bin_hi histo[i*2+1].wave = data[last_row].bin_lo histo[i*2].flux = data[n[i]].grp_data / (histo[i*2].wave - histo[i*2+1].wave) histo[i*2+1].flux = data[n[i]].grp_data/ (histo[i*2].wave - histo[i*2+1].wave) endfor return, histo end ;==================== pro custom_group, data, out=out, file_ps=file_ps ;GROUPING: Equal to 1 if the row is the first in a group; equal to -1 otherwise. ;QUALITY: The quality flag for a given group (same conventions as for the ftool GRPPHA): 0 for good (grouped) data; 5 for data labeled as bad by the user (within a tab), and 2 for data labeled as questionable by dmgroup (incomplete groups, etc.). ;GRP_NUM: The number of the group in which that row falls. ;CHANS_PER_GRP: The total number of rows in the given group. ;GRP_DATA: The total number of counts in the given group. ;GRP_STAT_ERR: The statistical error for the group. If the error column is supplied, it is the value of these errors added in quadrature. If no error column is supplied, it is sqrt(grp_data). data = mrdfits(data, 1, head) data_new = data ndata = n_elements(data) O7_plot = [20,23] Ne9_plot = [13,14] Mg11_plot = [9,9.6] Si13_plot = [6.4,7] bin = fltarr(2,4) bin[*,0] = [21.45,22.5] ; O 7 bin[*,1] = [13.38,13.75] ; Ne 9 bin[*,2] = [9.15,9.35] ; Mg 11 bin[*,3] = [6.6,6.8] ; Si 13 lim = fltarr(2,4) for i=0, 3 do begin group = where(data.grouping eq 1, ngroup) n = where( (data[group].bin_hi le bin[1,i]) and data[group].bin_hi gt bin[0,i], nc ) first_group = group[n[0]] first_row = group[n[0]] last_group = group[n[nc-1]] last_row = last_group + data_new[last_group].chans_per_grp n_in_new_group = last_row - first_row counts_in_new_group = total( data_new[group[n]].grp_data ) err_counts_in_new_group = sqrt(total( data_new[group[n]].grp_stat_err^2 )) lim[0,i] = data_new[last_row].bin_lo lim[1,i] = data_new[first_row].bin_hi data_new[group[n[1:nc-1]]].grouping = -1 data_new[first_row:last_row-1].grp_num = data_new[first_group].grp_num data_new[first_row:last_row-1].grp_data = counts_in_new_group data_new[first_row:last_row-1].chans_per_grp = n_in_new_group data_new[first_row:last_row-1].grp_stat_err = err_counts_in_new_group data_new[first_row:last_row-1].quality = 0 data_new[last_row:ndata-1].grp_num -= (nc-1) ;print, [transpose(data[first_row-1:last_row].grp_num), transpose(data[first_row-1:last_row].quality), transpose(data_new[first_row-1:last_row].quality) ] ;print ;print, data[group[n]].grp_num ;print, data[group[n]].chans_per_grp ;print, data[group[n]].grp_data, '=>', counts_in_new_group ;print, data[group[n]].grp_stat_err, '=>', err_counts_in_new_group ;print, first_row, last_row, '=>', n_in_new_group ;print, data[first_group].grp_num, data[last_group].grp_num, '=>', nc ;help, data[O7_first] ;help, data[O7_first+1] endfor if keyword_set(out) then mwrfits, data_new, out, head, /create data_histo = pha_grp_histo(data) data_new_histo = pha_grp_histo(data_new) if keyword_set(file_ps) then begin set_plot,'PS' device, /portrait, /encaps, /color, /inches, xsize=8, ysize=8, file=file_ps endif !p.multi=[0,2,2] load_color_vp plot_vp2 plot, data_histo.wave, data_histo.flux, xrange = O7_plot, title='O7' oplot, replicate(21.605, 2), !y.crange, line=1 oplot, replicate(21.804, 2), !y.crange, line=1 oplot, replicate(22.1, 2), !y.crange, line=1 ; ? oplot, replicate(lim[0,0], 2), !y.crange, line=1, color=2 oplot, replicate(lim[1,0], 2), !y.crange, line=1, color=2 oplot, data_new_histo.wave, data_new_histo.flux, color=3 plot, reverse(data_histo.wave), reverse(data_histo.flux), xrange = Ne9_plot, title='Ne9' oplot, replicate(13.447, 2), !y.crange, line=1 oplot, replicate(13.553, 2), !y.crange, line=1 oplot, replicate(13.699, 2), !y.crange, line=1 oplot, replicate(lim[0,1], 2), !y.crange, line=1, color=2 oplot, replicate(lim[1,1], 2), !y.crange, line=1, color=2 oplot, data_new_histo.wave, data_new_histo.flux, color=3 plot, reverse(data_histo.wave), reverse(data_histo.flux), xrange = Mg11_plot, title='Mg11' oplot, replicate(9.231, 2), !y.crange, line=1 oplot, replicate(9.314, 2), !y.crange, line=1 oplot, replicate(9.169, 2), !y.crange, line=1 oplot, replicate(lim[0,2], 2), !y.crange, line=1, color=2 oplot, replicate(lim[1,2], 2), !y.crange, line=1, color=2 oplot, data_new_histo.wave, data_new_histo.flux, color=3 plot, reverse(data_histo.wave), reverse(data_histo.flux), xrange = Si13_plot, title='Si13' oplot, replicate(6.647, 2), !y.crange, line=1 oplot, replicate(6.687, 2), !y.crange, line=1 oplot, replicate(6.740, 2), !y.crange, line=1 oplot, replicate(lim[0,3], 2), !y.crange, line=1, color=2 oplot, replicate(lim[1,3], 2), !y.crange, line=1, color=2 oplot, data_new_histo.wave, data_new_histo.flux, color=3 if keyword_set(file_ps) then device, /close return end