Subject: Re: grppha From: VĂ©ronique Petit Date: 6/15/12 1:36 PM To: David Cohen CC: "Gagne, Marc" , "Leutenegger, Maurice A. (GSFC-662.0)[OAK RIDGE ASSOCIATED UNIVERSITIES (ORAU)]" , Jacob Ellington Neely Hi David, What I did for tau sco and Ori2A is that I did the fancy (adaptative snr) grouping in ciao, then I modified the grouping directly in the pha file with IDL to group the He-like complex. I am attaching the IDL code. Main code is custom_group, plus some other codes I think I call, although it might call some other of my other codes, so tell me if something is missing.. Alternatively, you can do the SNR grouping you want and send me the pha file. Cheers, V. On 12-06-15 11:08, David Cohen wrote: > X-ray dudes, > > Jake and I are doing some apec modeling of zeta Ori and zeta Pup's RGS spectra. > > We're binning up the He-like complexes, since apec doesn't handle the altered f/i ratios. And we'd also like to bin the continuum to meet a minimum count or S/N criterion. > > It seems like grppha won't let us apply both types of grouping (group Nmin Nmax bins; and group min value) sequentially to a given pha file. > > Does that make sense to you? Do you think it's possible to apply both types of grouping criteria? > > Thanks in advance for your advice. > > David custom_group.pro 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 pha_grp_histo.pro function pha_grp_histo, data ndata = n_elements(data) ; nt = n_tags(data) ; if nt le 7 then begin ; for i=0, ngroup-1 do begin ; histo[i*2].wave = data[].bin_hi ; histo[i*2+1].wave = data[last_row].bin_lo ; endif else begin 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 ; endelse return, histo end read_hetg_xspec.pro function read_hetg_xspec, file nlines = file_lines(file) ;print, nlines data_parse = strarr(nlines) get_lun, lun openr, lun, file readf, lun, data_parse close, lun free_lun, lun ;print, transpose(data_parse) search = strpos(data_parse, '!') n = where(search ge 0, ncount) if ncount lt 0 then begin print, 'problem with file header...' return, -1 endif print, 'first line of data at', n[0]+2 data_parse = data_parse[n[0]+1:nlines-1] nlines -= (n[0]+1) ; how many columns in data search = strsplit(data_parse[0], count=ncol) print, data_parse[0] print, 'there are ', ncol, ' column' struct = {wave:0D, bin:0D, flux:0D, err_I:0D} if ncol gt 4 then begin name = 'model'+string(indgen(ncol-4)+1, format="(I0)") help, name help, replicate(0D, ncol-4) for i=1, (ncol-4) do struct = create_struct( struct, 'model'+string(i, format="(I0)"), 0D ) endif search = strpos(data_parse, 'NO') n = where(search ge 0, ncount) nspec = ncount+1 print, 'there are ', nspec, ' spectra' if nspec eq 1 then spec_pos = [-1, nlines] else spec_pos = [-1, n, nlines] print, spec_pos for i=0, nspec-1 do begin nwave = spec_pos[i+1]-1 - spec_pos[i] print, i, ' : ', nwave spec_data = replicate(struct, nwave) reads, data_parse[ spec_pos[i]+1 : spec_pos[i+1]-1 ], spec_data ;spec_data = reverse(spec_data) ; note the wave scale is reversed if i eq 0 then data = create_struct('spec1', spec_data) else data = create_struct(data, 'spec'+string(i+1, format="(I0)"), spec_data) ; for j= spec_pos[i]+1, spec_pos[i+1]-1 do begin endfor return, data end ; col 1 = wavelength ; col 2 = bin half-width ; col 3 = data ; col 4 = uncertainty ; col 5 = model Attachments: custom_group.pro 5.1 KB pha_grp_histo.pro 780 bytes read_hetg_xspec.pro 1.6 KB