Hi Raul, Do you have any zeros in your data? If you're calculating the fit statistic based on your data, and some have errors = 0, that can lead to large chi-squared. Try setting fstat_mod and refitting! Cheers, Joey Sent from my iPhone On May 24, 2013, at 3:40 PM, Raul Eduardo Puebla <rep54_at_email.domain.hidden > Hi Dave, > > > I tried to do the fit. I fixed the densities, and constrain the other > parameters as follows (data_vec2 corresponds to MEG +- 1 data): > > fit_fun( "xaped(1) * phabs(1)" ); > list_par; > set_par("xaped(1).norm",0.0005,0); > set_par("xaped(1).temperature",1e+7,0, 1.e6, 1.e8); > set_par("xaped(1).density",1,1); > use_thermal_profile; > set_par( "xaped(1).vturb", 300, 0, 0.0, 1000); > set_par( "phabs(1).nH", 0.1, 0, 0, 0.5); > list_par; > renorm_counts; > list_par; > exclude( all_data ); > list_data; > include(data_vec2); > list_data; > group_data( data_vec2, 4 ); > xnotice( data_vec2, lo, lf ); > fit_counts; > popt.res=4; > plot_counts( -data_vec2,popt ); > list_par; > > > But it get the following output for model parameters: > > Parameters[Variable] = 7[4] > Data bins = 1052 > Chi-square = 2.04803e+13 > Reduced chi-square = 1.954227e+10 > Graphics device/type (? to see list, default /xwin): > Invalid: plot limits > plot failed > xaped(1) * phabs(1) > idx param tie-to freeze value min > max > 1 xaped(1).norm 0 0 -5.087167e-28 0 > 0 > 2 xaped(1).temperature 0 0 1e+07 1000000 > 1e+08 > 3 xaped(1).density 0 1 1 0 > 0 > 4 xaped(1).vturb 0 0 300 0 > 1000 > 5 xaped(1).redshift 0 1 0 0 > 0 > 6 xaped(1).metal_abund 0 1 1 0 > 0 > 7 phabs(1).nH 0 0 0.1 0 > 0.5 > > > When I fix the xaped(1).norm to 0.0005, I get: > > Parameters[Variable] = 7[3] > Data bins = 1052 > Chi-square = 2.86434e+54 > Reduced chi-square = 2.730544e+51 > Graphics device/type (? to see list, default /xwin): > xaped(1) * phabs(1) > idx param tie-to freeze value min > max > 1 xaped(1).norm 0 1 0.0005 0 > 0 > 2 xaped(1).temperature 0 0 1000000 1000000 1e+08 > 3 xaped(1).density 0 1 1 0 > 0 > 4 xaped(1).vturb 0 0 0 0 > 1000 > 5 xaped(1).redshift 0 1 0 0 > 0 > 6 xaped(1).metal_abund 0 1 1 0 0 > 7 phabs(1).nH 0 0 0.5 0 > 0.5 10^22 > > > I the both cases, I get a model with counts=0, for all channels. I think > I am still missing something. > > > Thanks > > Raul > > > > > > > > On Wed, May 22, 2013 4:50 pm, David P. Huenemoerder wrote: >> >> Raul> I have combined the data from HEG +/- 1 and MEG +-1. Is it >> Raul> possible to fit the combined spectrum with some xpec model (or >> Raul> what else), using >> >> Raul> fit_fun("apec"); >> Raul> fit_counts; >> >> Hi Raul, >> >> Short answer is "yes". But there is a lot implied in your email. >> Here's what I typically do (or the equivalent of, and there are, of >> course, variations on the theme): >> >> Load your data, arfs, and rmfs: >> >> h = load_data( "pha2", [3,4,9,10] ); % load HEG and MEG rows >> >> amm = load_arf( "meg_-1.arf" ); >> amp = load_arf( "meg_1.arf" ); >> ahm = load_arf( "heg_-1.arf" ); >> ahp = load_arf( "heg_1.arf" ); >> >> rmm = load_arf( "meg_-1.rmf" ); >> rmp = load_arf( "meg_1.rmf" ); >> rhm = load_arf( "heg_-1.rmf" ); >> rhp = load_arf( "heg_1.rmf" ); >> >> If you use the tgcat.sl suite and files donwloaded from tgcat.mit.edu, >> all that can be done via: >> >> d = load_set_acis( "your_data_directory", [3,4,9,10] ); >> >> (it assumes the tgcat generic file naming scheme) >> >> (I like to keep the returned values, especially for use in scripts. For >> interactive use, you don't necessarily need to use "x = ...".) >> >> Set flags for combining the data: >> >> gh = combine_datasets( h[[0,1]] ); >> gm = combine_datasets( h[[2,3]] ); >> >> (without the variable h, you could do: gh = combine_datasets( [1,2] ); >> etc) >> >> Note that this does not add any spectra, arfs or rmfs, but solely sets >> flags identifying which are to be combined. >> >> You can use the xspec apec model if you want. If you have AtomDB >> installed, then there are advantages to using a native ISIS plasma >> model, since then you can query the database after model evaluation for >> the explicit line contributions (e.g., line ids, fluxes, etc, per >> wavelength range, by element, ion; that's another topic; see >> create_aped_fun, plot_group, for examples). >> >> For the question at hand, we'll use the xspec model: >> >> fit_fun( "apec(1) * phabs(1)" ); >> >> list_par; >> >> set_par( ...); >> >> exclude( all_data ); % in case I have multiple datasets loaded >> >> include( h ); % include the spectra of interest >> >> group_data( h, 2 ); % bin up a bit, uniformly (more or less if needed) >> >> xnotice( h, 1.7, 25 ); % notice the data and region of interest >> >> fit_counts; >> >> Fitting is the same whether you have done combine_datasets, or not. In >> this case, you have specified that MEG +-1 be combined, and that HEG +-1 >> be combined. When fit, the model is folded through the response for >> each grating and order, then before computing the statistic, the MEG >> data are added, the MEG model counts are added, and ditto for HEG, then >> the statistic computed (see the help for combine_datasets). So MEG and >> HEG are fit jointly, while each one is combined. If you are photon >> starved and want to combine them all, you first need to >> match_dataset_grids: >> >> match_dataset_grids( h[[2,0,1]] ); % put HEG on MEG grid >> >> g = combine_datasets( h ); >> >> >> Now, visualizing is another issue. If you use the tgcat.sl suite, it >> will load fancy_plots.sl, which has functions for plotting combined >> data: >> >> popt.res = 4; % combined residuals >> plot_counts( -[2,3], popt ); % negative sign means combine before >> plotting >> >> or equivalently: >> >> plot_counts( -gm, popt ); % using group id. >> >> >> There are also plot_data() (count rate) and plot_unfold() (fluxed >> units). >> >> >> >> >> -- Dave >> >> David Huenemoerder 617-253-4283 (o); -253-8084 (f); >> http://space.mit.edu/home/dph >> MIT Kavli Institute for Astrophysics and Space Research >> One Hampshire Street, Room NE80-6065 >> Cambridge, MA 02139 >> [Admin. Asst.: Elaine Tirrell, 617-253-7480, egt_at_email.domain.hidden>> >> > > ---- > You received this message because you are > subscribed to the isis-users list. > To unsubscribe, send a message to > isis-users-request_at_email.domain.hidden> with the first line of the message as: > unsubscribe > ---- You received this message because you are subscribed to the isis-users list. To unsubscribe, send a message to isis-users-request_at_email.domain.hiddenwith the first line of the message as: unsubscribeReceived on Fri May 24 2013 - 16:07:49 EDT
This archive was generated by hypermail 2.2.0 : Thu Jun 06 2013 - 12:24:50 EDT