; ;+ ; add a line to the subtitle of the spectral plots ; more_info must be set to display this ;- pro more_subtitle, subtitle, sigmaa=sigmaa, chisqr=chisqr, more_info=more_info if keyword_set(sigmaa) and keyword_set(chisqr) and keyword_set(more_info) then $ subtitle = subtitle+ $ '!c Red. Chisqr: '+string(chisqr,form='(g8.3)') + $ ' Uncertainties: '+ $ arr2str(delim=', ',strcompress(/remove, string(sigmaa,form='(g8.3)') )) end ;+ ; ; NAME: spec_plot ; ; PURPOSE: Plot an observed energy spectrum and its model ; ; CATEGORY: SPEX, GRAPHICS ; ; ; CALLING SEQUENCE: ; spec_plot, edges=edges, drm=drm, e_in=e_in, apar=apar, $ ; conversion=conversion, pconv=pconv, cflux=cflux, ecflux=ecflux, units=units,$ ; set=set, xtype=xtype, ytype=ytype, colors=colors, uncert=uncert, cal_test=cal_test ; CALLS TO: Fcount_rate ; ; ; KEYWORD INPUTS: ; edges: 2 x Num_channels energy edges ; drm: Detector Response defining relationship between input model and output response ; DRM is Num_channels x Num_photon ; e_in: Input energy edges for photon spectrum, 2 X Num_photon bins ; apar: Parameters of model ; sigmaa: Uncertainties on model parameters ; cflux: Counting rate flux, counts/uflux, often counts/cm2/s/keV, Num_channels long. ; ecflux: Uncertainties on cflux ; bflux: Background flux ; ebflux: Uncertainty on bflux ; units: Units of cflux, expressed as needed for figure, often' s!u-1!n cm!u-2!n keV!u-1!n' ; set: if set then interactively select fitting parameters ; interval: Interval number used to reference a series of fits. ; Trange: Time range covered by observation, for subtitle. ; Chisqr: Value of reduced chi-square for this fit. ; ; cal_test 1, scaled plot comparing net source count/rate and model count/rate ; In the cal_test mode we are trying to highlight the correspondance between ; the details of the predicted model spectrum and the observed spectrum looking ; for any obvious systematic problems in the gain/detector-response. ; Therefore, we plot the net source count rate spectrum over the range ; indicated by edges(*,wuse) and overplot with the predicted source count rate. ; Furthermore if a linear plot is desired then the plot is autoscaled by ; multiplying by a power-law which will best flatten the plot and the yrange ; is selected from extrema of the product of the observed and the power-law over the ; the channels in wuse. ; Title -Title used on figure, string. ; xtype -If set log xaxis, otherwise linear ; ytype -If set log yaxis, otherwise linear ; Colors - Vector of colors to use, referenced to Linecolors.pro ; Nominally [5,9,3,4,7] for ; linecolor 5 (yellow) source counts ; linecolor 9 (sky blue) folded model counts, ; linecolor 3 (scarlet) total continuous model, ; linecolor 4 (orange) background counts, ; linecolor 7 (green) continuous model components. ; group=group ; The new edges are the low edge of the first channel in each group and the low ; edge of the second channel. ; uncert - minimum fractional uncertainty on background subtracted flux ; pconv - pileup correction to conversion factor, optional, see run_curvefit ; OVERPLOT - If set, overplot new spectrum ; RESIDUALS - If set, plot residuals. ; OUTPUTS: ; conversion ; ; OPTIONAL OUTPUTS: ; none ; ; COMMON BLOCKS: ; Function_com ; ; SIDE EFFECTS: ; none ; ; RESTRICTIONS: ; none ; ; PROCEDURE: ; If 'COUNTS' appears in units, flux is displayed with no additional scaling, ; and the bflux is also plotted. Otherwise, the (cflux-bflux) is divided by the ; conversion factor (see GET_CONVERSION.PRO) computed from the ; drm and model function found in function_com. ; ; MODIFICATION HISTORY: ; ras, 1994 ; ras, 9-Mar-1995, only 6 parameters per line in the subtitle ; ras, 22-mar-95, added cal_test keyword ; ras, 5-apr-95, fixed bug in energy bin averaging, creating gmatrix ; rasm 24-jan-96, added residuals plot ;- pro spec_plot, edges=edges, delta_light=delta_light, drm=drm, e_in=e_in, apar=apar, group=group, $ interval=ix, cflux=cflux, ecflux=ecflux, bflux=bflux, ebflux=ebflux, trange=trange, $ set=set, units=units, xtype=xtype, ytype=ytype, pconv=pconv, $ wuse=wuse, yrange = yrange, xrange = xrange, cal_test=cal_test, residuals=residuals, $ sigmaa=sigmaa, chisqr=chisqr, more_info=more_info, $ uncert=uncert, colors=colors, title=title, overplot=overplot @function_com pmulti = !p.multi checkvar, residuals, 0 if residuals then !p.multi=[0,1,2] checkvar, colors, [5,9,3,4,7] if n_elements(colors) eq 1 then colors = [colors,9] if n_elements(colors) eq 2 then colors = [colors,3] if n_elements(colors) eq 3 then colors = [colors,4] if n_elements(colors) eq 4 then colors = [colors,7] ;If fitting parameters are passed, then show the model function fit = keyword_set(total(abs(fcheck(apar,0.0)))) ;if in the units passed it says COUNTS, then plot in counts ;otherwise show a photon y axis, and project data into photons. checkvar,units,'Counts s!u-1!n cm!u-2!n keV!u-1!n' if strpos(strupcase(units),'COUNT') ne -1 then count=1 else count=0 cal_test = keyword_set(cal_test) if cal_test then begin ;fit = 0 ;won't show photon model directly ;count = 1 ;scaled count units endif ;Do you want to set the parameters interactively? set = fcheck(set, 0) ;log-log is default plot style checkvar, xtype, 1 checkvar, ytype, 1 spxrange = fcheck(xrange, [0.0,0.0]) spyrange = fcheck(yrange, [0.0,0.0]) checkvar, uncert, 0.0 nchan = n_elements(cflux) ;number of ungrouped channels wuse = fcheck(wuse, indgen(nchan)) ecflux1= fcheck( ecflux, cflux*0.0) bflux1 = fcheck( bflux, cflux*0.0) * count ;no background shown in photon plot ebflux1 = fcheck( ebflux, cflux*0.0) * count ;no background shown in photon plot cflux1 = cflux - fcheck( bflux, cflux*0.0) ;cflux is assumed background subtracted ecflux1 = sqrt( ecflux1^2 + fcheck(ebflux,cflux*0.0)^2) if keyword_set(trange) then begin ix = fcheck(ix,0) sapar = strcompress(/remove, string(apar,form='(g9.4)')) numpar = n_elements(apar) apar_list = arr2str(delim=', ',sapar(0:5<(numpar-1)) ) if numpar gt 6 then for ii=6,numpar-1,6 do $ apar_list = apar_list+'!c'+arr2str(delim=', ',sapar(ii:(ii+5)<(numpar-1)) ) subtitle = 'Interval '+strtrim(ix,2)+'!c' + $ strmid(atime(trange(0),/pub),9,12) +' - '+$ strmid(atime(trange(1),/pub),9,12) +$ '!c' +f_model+' parameters: '+ apar_list endif more_subtitle, subtitle, sigmaa=sigmaa, chisqr=chisqr, more_info=more_info edge_products, edges, wid=wedges, mean=emdge ;------------------------------------------------- ;This code bins the channels together as specified by the variable group, if it exists ;The new edges are the low edge of the first channel in each group and the low ;edge of the second channel if keyword_set(group) then begin nchan1 = n_elements(group(0,*)) ;number of grouped channels ;If the grouping extends to the last channel, then the edge array needs to be extended gedges = reform([edges(*),edges(*,nchan-1)],2,nchan+1) edges1 = transpose( reform( [ (gedges(0,group(0,*)))(*), (gedges(0,group(1,*)))(*)],nchan1,2)) edge_products, edges1, wid= wedges1 ;18-aug-94, ras, allow for difference in delta bin edges and bin width used for count_2_flux ;bin edges given by edges ;normalization width from pulse-height amplitudes given by delta_light ; wdelta = wedges1*0.0 ; for i=0,n_elements(wdelta) do wdelta(i) = $ ; total(delta_light(group(0,i):group(1,i)-1)>group(0,1)) gmatrix = fltarr(nchan1,nchan) for i=0,nchan1-1 do gmatrix(i,group(0,i):group(1,i)-1>group(0,1)) = $ f_div(delta_light(group(0,i):group(1,i)-1>group(0,1) ), $ total(delta_light(group(0,i):group(1,i)-1>group(0,1)))) erange = limits( edges(*,wuse)) wuse1 = where( edges1(1,*) gt erange(0) and edges1(0,*) lt erange(1)) if fit then drm1 = gmatrix# drm cflux1 = gmatrix# cflux1 bflux1 = gmatrix# bflux1 ecflux1 = sqrt(gmatrix^2 # ecflux1^2) endif else begin edges1 = edges wedges1 = wedges if fit then drm1 = drm wuse1 = wuse endelse edge_products, edges1, mean=emdge1 if fit then begin edge_products, e_in, mean = em, width=w_in ord_em = sort(em) get_conversion, e_in=e_in, w_in=w_in, drm=drm1, e_out=edges1, $ f_model=f_model, apar=apar, conversion=conv, model_in=model, model_out=model_out conv = f_div(conv, fcheck(pconv,1.0)) crate = model_out * conv ; crate = fcount_rate(em=em,model=model, drm=drm1, e_in=e_in, par=apar) ; conv = crate(*) / (call_function(f_model, emdge1, apar))(*) mflux = crate(*) endif conv = fcheck(conv, cflux1*0.+1.) if count then conv = fltarr(n_elements(cflux1))*0.0 + 1. if cal_test then begin countsmod_plot, edges, e_in, drm, cflux, bflux, xrange=spxrange, apar=apar, $ colors=colors, yrange=spyrange, f_model=f_model, wuse=wuse, delta_light=delta_light goto, leavehere endif flux = cflux1/conv eflux = (ecflux1 > uncert*cflux1)/conv if spyrange(0) eq 0 and spyrange(1) eq 0 then $ spyrange = limits( f_div((flux+bflux1)(wuse1),conv(wuse1))) > 1e-6 < 1e4 if spxrange(0) eq 0 and spxrange(1) eq 0 then $ spxrange = limits( edges1(*,wuse1)) ;insert block for testing gain uncertainties common gain_mod_test, gmod, win11, win12 checkvar, gmod, 1.0 if not keyword_set(overplot) then begin plot, edges1(0,*)*gmod, (flux+bflux1)/gmod, psym=3, xtype=xtype, ytype=ytype, ytitl=units, $ xtitle = 'Energy (keV)', subtitle=fcheck(subtitle,''),$ xrange=spxrange, yrange=spyrange,title=title ;Save for future overplotting win11 = {winsav, x:!x, y:!y, clip:!p.clip} endif else begin winold = win11 ;set the axes x=winold.x y=winold.y position = [winold.x.window(0),winold.y.window(0), $ winold.x.window(1),winold.y.window(1)] if y.type eq 1 then ycrange=10^y.crange else ycrange=y.crange if x.type eq 1 then xcrange=10^x.crange else xcrange=x.crange plot,xcrange,ycrange,/nodata,xstyle=5,ystyle=21,/noerase, $ position=position,ytype=y.type,xtype=!x.type !p.clip=winold.clip endelse datplot, 0,0, xs=edges1*gmod, (flux+bflux1)/gmod, color=fcolor(colors(0)) if total(abs(bflux1)) ne 0.0 then $ ;plot the background in count or flux units datplot, 0,0, xs=edges1*gmod, bflux1/gmod, color=fcolor(colors(3)) if total(eflux) ne 0 then $ eplot, emdge1*gmod, (flux+bflux1)/gmod, ey=eflux, color=fcolor(colors(0)) if fit then begin datplot, 0,0, xs=edges1,f_div(mflux+bflux1,conv),color=fcolor(colors(1)) if keyword_set(set) * (not count) then begin widget=f_use_widget(/test,/continue ) printx,string(/print, 'Parameters for '+f_model+':',apar) buttons=['CHANGE PARAMETERS' , 'CONTINUE'] if widget then $ ans= respond_widg( title='SPECTRAL PLOT DIALOG BOX', $ message = [' The Photon Energy Spectrum for the Current Interval',$ ' will be plotted along with the model spectrum.', $ ' To change the parameters of the model, '+f_model + ',', $ 'Click on CHANGE PARAMETERS, otherwise CONTINUE'], /column,$ xoffset = FCHECK(xoffset,0), yoffset=FCHECK(yoffsets), $ buttons= buttons) $ else $ ans = emenu(buttons, $ mtitle= [' The Photon Energy Spectrum for the Current Interval',$ ' will be plotted along with the model spectrum.', $ ' To change the parameters of the model, '+f_model + ',', $ 'Select CHANGE PARAMETERS below, otherwise CONTINUE below.']) graphics_page if ans(0) eq 0 then begin;obtain new model parameters if widget then option_changer, 'Fit Parameters for '+f_model,'apar',direct=apar endif alpha_page oplot, em(ord_em), (call_function(f_model, e_in, apar))(ord_em), color=fcolor(colors(2)) endif else oplot, em(ord_em), model(ord_em), color = fcolor(colors(2)) ;Plot model for each individual component norm_param = model_components(f_model) if n_elements(norm_param) gt 1 then $ for i=0,n_elements(norm_param)-1 do begin apar_norm = apar apar_norm(norm_param) = 0.0 apar_norm(norm_param(i)) = apar(norm_param(i)) oplot, em(ord_em), (call_function(f_model, e_in, apar_norm))(ord_em), color=fcolor(colors(4)) endfor endif if not keyword_set(overplot) then timestamp if residuals then begin resid = f_div( flux-f_div(mflux,conv), eflux) res_range = minmax(resid(wuse1)) if not keyword_set(overplot) then begin plot, edges1(0,*)*gmod, resid, psym=3, $ xtype=xtype, ytitl='Normalized Residuals', $ xtitle = 'Energy (keV)',$ xrange=spxrange, yrange=res_range win12 = {winsav, x:!x, y:!y, clip:!p.clip} endif else begin winold = win12 ;set the axes x=winold.x y=winold.y position = [winold.x.window(0),winold.y.window(0), $ winold.x.window(1),winold.y.window(1)] if y.type eq 1 then ycrange=10^y.crange else ycrange=y.crange if x.type eq 1 then xcrange=10^x.crange else xcrange=x.crange plot,xcrange,ycrange,/nodata,xstyle=5,ystyle=21,/noerase, $ position=position,ytype=y.type,xtype=!x.type !p.clip=winold.clip endelse datplot, 0,0, xs=edges1*gmod, resid, color=fcolor(colors(0)) endif leavehere: !p.multi = pmulti alpha_page end