OSPEX Command Line Examples

Plotting the data spectrum for an interval in photon space after fitting:

emid = o->getaxis(/ct_energy, /mean)
ph_data = o->calc_summ(item='data_photon_flux', this_interval=2)
plot, emid, ph_data, /ylog,/xlog,yra=[1.e-2,1.e5],psym=10, /xstyle

Plotting the time profile for a single raw energy band:

d=o->getdata(spex_units='flux') ; could be 'counts', 'rate', or 'flux'
ut = o->getaxis(/ut, /mean)
utplot, anytim(ut,/ext), d.data[10,*],/ylog ; this plots energy bin number 10

Grouping energy or time bins:

d will contain structure of data from file in 'raw' (from file) energy and time bins. db will contain the grouped data in an array. To retrieve the binned errors and/or livetimes as well, use keywords.

d=o->getdata(spex_units='flux')   ; could be 'counts', 'rate', or 'flux'
dobj = o->get(/obj,class='spex_data')

eband = get_edges([3.,6.,12.,25.,50.], /edges_2)   ; or     eband = o->get(/spex_eband)   gets energy bands used in GUI
db = dobj->bin_data(data=d, intervals=eband)
or
db = dobj->bin_data(data=d, intervals=eband, eresult=eresult, ltime=ltime) ; to get grouped errors and livetime as well as data

tband = o->get(/spex_tband) ; time intervals used in GUI for grouping, or set your own time intervals [2,n]
db = dobj->bin_data(data=d, intervals=tband, /do_time_bin)    ; grouping in time intervals, NOTE the /do_time_bin keyword
or
db = dobj->bin_data(data=d, intervals=tband, /do_time_bin, eresult=eresult, ltime=ltime)

Grouping energy bins and plotting time profiles:

Below, d will contain structure of data from file in 'raw' (from file) energy and time bins. db will be an array containing d.data grouped into your new energy bins, dimensioned [n_eband,ntime]

eband = [[10.,20.],[20,30],[30,50]]
d=o->getdata(spex_units='rate') ; could be 'counts', 'rate', or 'flux'
dobj = o->get(/obj,class='spex_data')
db = dobj->bin_data(data=d, intervals=eband)

ut = o->getaxis(/ut, /mean)
utplot, anytim(ut,/ext), db[0,*], /ylog ; this plots grouped energy band 0
linecolors
outplot, anytim(ut,/ext), db[1,*], col=12
outplot, anytim(ut,/ext), db[2,*],col=7

Retrieving background-subtracted data, observed data, and background data grouped into analysis time intervals:

Assuming we have 77 raw energy bins, and 3 analysis time intervals:

d = o->getdata(class='spex_fitint', spex_units='flux') ; could be 'counts', 'rate', or 'flux'
help,d

** Structure <1ae44fa0>, 8 tags, length=7392, data length=7392, refs=2:
DATA FLOAT Array[77, 3]
EDATA FLOAT Array[77, 3]
LTIME FLOAT Array[77, 3]
OBSDATA FLOAT Array[77, 3]
EOBSDATA FLOAT Array[77, 3]
BKDATA FLOAT Array[77, 3]
EBKDATA FLOAT Array[77, 3]
BKLTIME FLOAT Array[77, 3]

Using plot_spectrum method:

Without PLOTMAN:

linecolors
o->plot_spectrum,this_interval=z,/show_fit,/use_fitted,spex_units='flux',/bksub,/overlay_back,$
xrange=[4.,1.e3],xstyle=1,/photons,/no_plotman,dim1_color=[0,2],fit_color=[6,3,5],charsize=1.4

With PLOTMAN, no legends, your own title:

o->plot_spectrum,/bk_sub,/overlay_back,/show_fit,/show_err,/use_fitted, spex_units='flux', fit_legend_loc=0,id='Count Flux, 19-January-2005 flare', legend_loc=0

Plotting items from the fit results:

o->plot_summ,xplot='time',yplot='flux_at_e',func_energy=[25.,50.]

Finding and plotting the Fermi GBM B0 data for a time interval in 2 energy bands:

time_range = ['23-Oct-2012 02:54', '23-Oct-2012 03:42']
gbm_find_data,date=time_range, pattern='cspec', det='b0', /copy, file=file
obj = ospex(spex_specfile=file,/no)
obj->set,spex_accum_time = time_range
gbmstruct = obj->getdata(spex_units='flux')
eband = [[200.,1000],[1000.,2000.]]
gbm = obj->bin_data(data=gbmstruct,intervals=eband,units_str=obj->getunits(class='spex_data'))
tgbm = obj->getaxis(/ut,/mean)
linecolors
utplot, anytim(tgbm,/ext), gbm[0,*], yra=[.005,.1],/ylog
outplot, anytim(tgbm,/ext), gbm[1,*], col=2

Using the spex_data class directly:

(Note - you can extract the spex_data object from an existing ospex object via dobj = o->get(/obj,class='spex_data') or create a standalone spex_data object as below.)

dobj = obj_new('spex_data')
o->set, spex_specfile= 'C:\Analysis\working\thick2_norm2\hsi_spectrum_d4.fits'
o->plot

o->plot, /pl_energy
o->plot, /pl_time

o->plot, /pl_time, intervals=[[3,6],[6,12],[12,25]]
o->plot, /pl_time, intervals=[[3,6],[6,12],[12,25]], dim1_sum=0 ; to plot the channels separately

o->plotman
o->plotman, /pl_time, intervals=[[3,6],[6,12],[12,25]], dim1_sum=0

d=o->getdata()
z=o->bin_data(data=d, intervals=[[3,6],[6,12],[12,25]]) ; group into energy bins
plot,z[0,*],/ylog,yra=[1.,1.e6],psym=10
oplot,z[1,*]
oplot,z[2,*]

Iain Hannah has a comprehensive example on GitHub that:

restores NuSTAR data from a save file
inserts the data into OSPEX using the 'SPEX_USER_DATA' method
fits the data, and
plots the count rate and count flux overplotted by the model.

Run the monte_carlo method and the chi2_map methods and plot the results (used to determine the uncertainties in parameter values):

o -> monte_carlo, niter=10, interval=0,savefile='test_mc.sav'
restore,file='test_mc.sav'
spex_monte_carlo_results, savefile='test_mc.sav', /plot, out_struct=mcs

o -> chi2_map, param_index=0, ntrials=100,savefile='test_c2m0.sav'
spex_chi2_map_results,savefile='test_c2m0.sav',/plot,out_struct=c2m_str
o -> chi2_map, param_index=1, ntrials=100,savefile='test_c2m1.sav'
spex_chi2_map_results,savefile='test_c2m1.sav',/plot,out_struct=c2m_str

Kim Tolbert
Last Modified: 10 November 2015