Hydra A - Event-2D Demo
Back to Event-2D Analysis page
dd@space.mit.edu
Hydra A - Event-2D DemoBack to Event-2D Analysis page |
| 11/13,29/07 dd@space.mit.edu |
The command window text is given below along with any plots or images created interspersed
appropriately.
The user-typed input is shown in bold and commentary is in this color.
hydra> .source ha_zodata Reading file: obs_4969_evt1_s3.fits - - - Reading file: obs_4969_evt1_s3.fits - - - Reading file: obs_4970_evt1_s3.fits - - - Reading file: obs_4970_evt1_s3.fits - - - - Data: e2d_meta[id].; e2d_data[id].; e2d_view[id]. id Mission-grat-instrument Date 2D-axes ignore prt m src events(k) exp(ks) target 10 CHANDRA-NONE-ACIS 2004.425 X,Y 0 - - - 1486.48 195.7 Hydra A 11 CHANDRA-NONE-ACIS 2004.425 RADIUS,ENERGY 0 - - - 1486.48 195.7 Hydra ALook at the datasets...
hydra> e2d_view_data(10);
e2d_view_data ( [ -150 , 150 , 2 ] , [ -150 , 150 , 2 ] );
total events displayed: 1.31773e+06 e2d_view[e2di].scale = 200%
bins w/events have 50%, 90%, and max numbers of: 30, 112, 1691
- - -
hydra> e2d_view_data(11);
e2d_view_data ( [ -5 , 150 , 2 ] , [ 0.1 , 8 , 0.075 ] );
total events displayed: 1.23637e+06 e2d_view[e2di].scale = 283%
bins w/events have 50%, 90%, and max numbers of: 38, 516, 2061
 |
For
dataset 11 the axes are RADIUS (left-right is -5 to 150 arc sec.) and ENERGY (vertically from 0.1 to 8 keV.) We won't explore this RADIUS-ENERGY space in the following, but it is certainly a useful way to compare and fit data and model for cluster's with a spectral variation primarily in the radial direction. |
- - -
hydra> .source ha_model_simple
Spectrum created from: ha_3p5keV0p5Abund.par
- - -
Spectrum created from: ha_AGN.par
- - -
o check/update foreground NH ...
o check/update spectra ...
o check/update the component geometry ...
The model-defining file, ha_model_simple.sl, ends with an
s3d_list_model command which produces the following output to
summarize the model information and the "s3d_ps" user-defined model
parameters:
--- Source-3D Model ---
Model name : Hydra A - simple
Distance : 216000 kpc. Date : 2004.400
Foreground NH : 0.0494 x10^22 /cm^2
Cube radius : 100 arc sec with 47 cells/radius, cell size of 2.12766 arc sec.
- - -
s3d_NHarf_method = 0 and 'slope = 3
Nominal spectral grid : 1 A to 40 A, with delta-log of 0.0003
- - -
- Spectra: s3d_spec[id].
id min - max es ph/ks/100cm^2 file
1 0.31 - 12.40 64.7776 ha_3p5keV0p5Abund.par
19 0.31 - 12.40 1838.66 ha_AGN.par
- - -
- Components: s3d_comp[id].
id name ignore type norm ispec vtype velval vturb rndint rndclr
1 Beta w/cavity 0 xray 22.7 1 -1 0 0 500 [1,0,0]
19 AGN 0 xray 0.00316 19 -1 0 0 1 [0,1,1]
- - -
The s3d_ps user parameter structure values:
{rAGN=1, rbeta=13.5, beta=0.433, ratio=1, betax=-3.3, betay=0}
- - -
hydra> s3d_update;
o check/update foreground NH ...
o check/update spectra ...
o check/update the component geometry ...
Look at the model geometry in both a solid-3D view
and as the projection of components. For the simple model here these are
not very exciting looking - so here are the commands without the images -
see the E0102 demo for better examples.hydra> s3d_view_comps; hydra> s3d_proj_comps;Show the spectra and NH curve...
hydra> s3d_plot_spectrum(19); % 19 is the AGN spectrum Graphics device/type (? to see list, default /XSERVE): CR hydra> s3d_oplot_spectrum(19,4); % make it blue hydra> s3d_oplot_spectrum(1,2); % cluster emission in red- - - hydra> s3d_plot_nh;
- - -
This next routine, e2d_pars2fom, does the whole parameters-to-figure-of-merit process - essentially everything inside a "function evaluation" that a fitting routine would access.
The routine s3d_xray_photons generates Monte-Carlo photons from the current s3d model components and their assigned spectra.
The routine e2d_inst_xray takes the sky photons and creates detector
photons - a very poor man's MARX.
In future this will also do XMM, Suzaku, etc.
hydra> e2d_pars2fom;
Event2D: parameters to FOM
- update parameters [e2d_update_pars]
- updating the model ... [s3d_update]
o check/update foreground NH ...
o check/update spectra ...
o check/update the component geometry ...
- Dataset 10: generating photons from the model... [s3d_xray_photons]
o Beta w/cavity: generating 438119 MC photons...
o AGN: generating 1536 MC photons...
o Total generated = 439655
10: pass photons through instrument simulation... [e2d_inst_xray]
10: load simulation into model arrays [e2d_inst2mod]
- Dataset 11: generating photons from the model... [s3d_xray_photons]
o Beta w/cavity: generating 438119 MC photons...
o AGN: generating 1536 MC photons...
o Total generated = 439655
11: pass photons through instrument simulation... [e2d_inst_xray]
11: load simulation into model arrays [e2d_inst2mod]
- calculate and show the figure of merit... [e2d_list_fom]
o Data-Model 10, chi^2 = 215444 from 22448 cells; chi^2/cell = 9.59748; sum(D)/sum(M) = 1.50058
o Data-Model 11, chi^2 = 157788 from 7659 cells; chi^2/cell = 20.6016; sum(D)/sum(M) = 1.40783
o Total FOM = 373232
- - -
Look at the Data, Model and Residual images...
hydra> e2d_view_resid(10);
blue indicates: chi < -2 [ model is much greater ]
green indicates: -2 < chi < 0 [ model is greater]
yellow-green indicates: 0 < chi < 2 [ data is greater ]
red indicates: chi > 2 [ data is much greater ]
At Left: Color-coded residual image - see print-out above for color coding
info.
At Right: The observed data, +/-150" range.
|  |
|
- - -Re-adjust the viewing (and FOM calculating) ranges to extend to only 100" to be within the simulation size.
hydra> e2d_set_view(10, [-100.0,100.0,1.0],[-100.0,100.0,1.0]); hydra> e2d_set_view(11, [-5.0,100.0,1.0],[0.1,8.0,0.075]);Look at the 1-D projections of the data along each axis of the two data sets, to get histograms vs X, Y, RADIUS, and ENERGY. Output plots not shown - see them for yourself !
hydra> e2d_plot_ks(10,1); hydra> e2d_plot_ks(10,2); hydra> e2d_plot_ks(11,1); hydra> e2d_plot_ks(11,2);The plots above indicated that the number of simulated events (the norm) is too low - we'll change that, but first... Make a K-S plot and evaluate the p-value for the data and model RADIUS distributions, dataset 11, axis 1. This statistic is not sensitive to the agreement in number of events - but even so this K-S test does tell us that the model and data shapes do _not_ agree to within statistics at the "no-way-at-all" level (p-value=0.)
hydra> e2d_plot_ks(11,1);
o Data-Model 11, K-S on RADIUS: Max |diff| is 0.0145491 at 30.1942 arcsec
ks_test2 statistic of 0.0145491 gives p-value of 0
- - -
hydra> e2d_ignore(11);
hydra> e2d_list_fom;
o Data-Model 10, chi^2 = 48336.6 from 39947 cells; chi^2/cell = 1.21002; sum(D)/sum(M) = 1.17811
o Total FOM = 48336.6
- - -
hydra> norm_adjust = e2d_norm_scale(1.0, 1.8, 0.02499);
- Minimum FOM of 40323 at norm factor = 1.19992
and calculate this adjustment more accurately
with:hydra> norm_adjust = e2d_norm_optimum(1.e-4); - Minimum FOM of 40316.4 at norm factor = 1.20653 - - -Adjust just the "Beta" component's norm by this factor (the AGN norm was set to agree with the AGN counts so leave it alone) and then re-evaluate the model with this new norm to check the FOM:
hydra> s3d_comp[1].norm;
22.7
hydra> s3d_comp[1].norm = s3d_comp[1].norm * norm_adjust;
hydra> s3d_comp[1].norm;
27.3882
hydra> e2d_pars2fom;
...
o Data-Model 10, chi^2 = 39693 from 39973 cells; chi^2/cell = 0.992995; sum(D)/sum(M) = 0.976795
o Total FOM = 39693 OK, that's better !
- - -
First, as a preliminary let's get a sense of how sensitive the FOM is to shifting the center of the Beta model in X and Y and in changing the Beta parameter. We'll setup to scan the parameter s3d_ps.betax and then also the parameter s3d_ps.beta by using the user-redefinable routine e2d_update_pars and the public scanning variables, e2d_scan_vals and e2d_iscan. The scanning range is set by human guess and we'll look at the plots that come out...
hydra> public define e2d_update_pars ()
{
s3d_ps.betax = e2d_scan_vals[e2d_iscan];
}
hydra> e2d_scan_par([ [-7.0:7.0:0.999] ], 1);
====================
Event2D: parameters to FOM
...
- FOM = 60020.3
- Scan parameter value was: 6.986
====== ======
====================
Minimum FOM is at scan value = -3.004 giving a FOM of 39310.5 +/- 235.705
with a best norm adjust of 0.989559
*** Parameter(s) set to this best-FOM value, BUT model not re-evaluated. ***
Fitting gives minimum at scan value = -2.51876 +/- 0.280121, with a FOM of 39611.1
====================
The plot below shows the scan results of
FOM vs betax for the scan values (black). The The green curve is a quadratic fit to the
measured points and it's minimum is at -2.52.
The betax minima reported above (the fit value and -3.004 which is the
closest of the scan values) are close to
the manual starting value of -3.3. A change in betax of 4 or 5 arc seconds
produces a clear change in FOM.
Now, set e2d_update_pars to scan the beta value (after putting the betax value back to manual value):- - -
hydra> s3d_ps.betax = -3.3;
hydra> public define e2d_update_pars ()
{
s3d_ps.beta = e2d_scan_vals[e2d_iscan];
}
hydra> e2d_scan_par([ [0.4:0.9:0.0499] ], 1);
====================
Event2D: parameters to FOM
...
- FOM = 271994
- Scan parameter value was: 0.899
====== ======
====================
Minimum FOM is at scan value = 0.5996 giving a FOM of 39969.9 +/- 2851.24
with a best norm adjust of 0.991001
*** Parameter(s) set to this best-FOM value, BUT model not re-evaluated. ***
Fitting gives minimum at scan value = 0.570917 +/- 0.0118932, with a FOM of 40142.5
====================
This plot shows a clear change in FOM with the beta value - the minimum is
very near the manually chosen value of 0.6, and a change of 0.2 in beta is
pretty significant.
hydra> e2d_cg_fit; % outputs useful summary help to setup a fit.
hydra> cg_p0 = [0.0, 0.0, 0.5];
hydra> cg_smand = [4.0, 4.0, 0.2];
hydra> public define e2d_update_pars ()
{
variable delta_p0;
delta_p0 = cg_smand * cg_dir * e2d_scan_vals[e2d_iscan];
s3d_ps.betax = cg_p0[0] + delta_p0[0];
s3d_ps.betay = cg_p0[1] + delta_p0[1];
s3d_ps.beta = cg_p0[2] + delta_p0[2];
}
So this will start with the beta model centered at X,Y = 0,0
and with a beta value of 0.5 and then iterate on these to minimize
the FOM...
hydra> e2d_cg_fit(1); % note that the "1" is a dummy argument
===========================
Starting CG iteration 1
===========================
...
[ Takes ~ 90 minutes running as a single thread on a
"Two AMD Opteron 2220 64Bit dual-core processors" machine ]
...
*** The criteria are satisfied. ***
All done with these iterations.
The CG fit creates an HTML file logging it's progress and including links
to the various 1D scans that are part of the fitting process so that the
user can monitor (by looking at HTML during the fitting) and evaluate
after the fact how the fitting process went. The output for
the fit above is linked here:
cg_log_Mon_Nov_12_14_21_22_2007.html
Well, it stopped iterating after the 5th iteration because
the criteria to stop were met: a small change in the new parameter vector value
from the previous value AND a change in FOM less than 0.7 times the
statistical FOM noise level. The result gives a
hydra> print(cg_p0); % confirm the cg_p0 starting values -1.44003 1.28023 0.573406 hydra> cg_smand = [3.0,3.0,0.01]; % reduce, especially the beta smand value hydra> e2d_cg_fit(1); % flip - the - switchOK, now after ~ 80 minutes we get:
I may have gone overboard in the other direction in this case
by reducing the beta smand value to 0.01. Alright, do it one
more time starting again at -1.44,1.28,0.573 and using [3.0,3.0,0.04]
for the smand values... ~ 55 minutes later we get:
cg_log_Mon_Nov_12_18_35_33_2007.html
which has well-matched derivatives and gives best fit values:
hydra> cg_dir=[0.,0.,0.]; % will evaluate at p0 hydra> e2d_pars2fom; hydra> e2d_view_resid(10);
Much of how the fitting process is carried out by e2d_cg_fit is based on the reality that the FOM values have variation even when the same parameters (and data) are input. This also sets a limit of how well the process can be expected to determine a minimum and hence how many iterations are needed before deciding to stop. The fit example given here has a very high signal-to-noise and these issues are not so apparent. In the case of fitting a more subtle feature, e.g., a cavity in the cluster, the noise in the Monte-Carlo approach will be more noticable. Let's see... :-)
hydra> .source ha_model_cavity hydra> s3d_ps.cv1rad=20.0; % cavity 1 radius hydra> s3d_ps.cv1xyz=[15.0, -25.0, 0.0]; % cavity 1 location in 3D hydra> s3d_ps.cv2rad=-1;s3d_ps.cv3rad=-1; % "turn off" other 5 cavities. hydra> s3d_ps.cv4rad=-1;s3d_ps.cv5rad=-1;s3d_ps.cv6rad=-1; hydra> s3d_update; hydra> s3d_comp[11].rndint=20.0; % adjust for pretty picture. hydra> s3d_proj_comps;Read in the data, use only the X-Y imaging data. Select a region of the data and binning to use that focusses on the feature we're fitting. Finally, adjust the component 1 norm for best agreement with the starting model.
- - -
hydra> .source ha_zodata
hydra> e2d_ignore(11);
hydra> e2d_view_data (10, [ -20 , 80 , 2 ] , [ -80 , 10 , 2 ] );
e2d_view_data ( [ -20 , 80 , 2 ] , [ -80 , 10 , 2 ] );
total events displayed: 405017 e2d_view[e2di].scale = 588%
bins w/events have 50%, 90%, and max numbers of: 80, 470, 1691
hydra> e2d_pars2fom;
hydra> s3d_comp[1].norm=s3d_comp[1].norm * e2d_norm_optimum(1.e-4);
- Minimum FOM of 8725.44 at norm factor = 1.08761
- - -
hydra> e2d_pars2fom;
...
o Data-Model 10, chi^2 = 9272.23 from 2250 cells; chi^2/cell = 4.12099; sum(D)/sum(M) = 0.969905
o Total FOM = 9272.23
- - -
hydra> e2d_view_resid(10);
blue indicates: chi < -2 [ model is much greater ]
green indicates: -2 < chi < 0 [ model is greater]
yellow-green indicates: 0 < chi < 2 [ data is greater ]
red indicates: chi > 2 [ data is much greater ]
One cavity in the southwest: starting residuals, data, and model. FOM =
9272
|  |
|
- - -Now set up the 3 fitting paramters (the cavity x,y and radius) with pretty large scale sizes, [20.0,20.0,20.0] and do the fit...
hydra> cg_p0 = Double_Type[3];
hydra> cg_smand = Double_Type[3];
hydra> cg_p0 = [15.0, -25.0, 20.0];
hydra> cg_smand = [10.0, 10.0, 3.0];
hydra> public define e2d_update_pars ()
{
variable delta_p0;
delta_p0 = cg_smand * cg_dir * e2d_scan_vals[e2d_iscan];
s3d_ps.cv1xyz[0] = cg_p0[0] + delta_p0[0];
s3d_ps.cv1xyz[1] = cg_p0[1] + delta_p0[1];
s3d_ps.cv1rad = cg_p0[2] + delta_p0[2];
}
hydra> e2d_cg_fit(1);
Here's the log of what the fitting did for about 3 hours:
hydra> cg_dir=[0.,0.,0.];
hydra> e2d_update_pars;
hydra> e2d_update_pars_NULL;
hydra> s3d_ps.cv1xyz[0]; % confirm the current values for cavity 1
40.3781
hydra> s3d_ps.cv1xyz[1];
-54.4197
hydra> s3d_ps.cv1xyz[2];
0
hydra> s3d_ps.cv1rad;
33.1628
hydra> e2d_pars2fom;
...
o Data-Model 10, chi^2 = 5134.72 from 2250 cells; chi^2/cell = 2.2821; sum(D)/sum(M) = 0.955752
o Total FOM = 5134.72
hydra> s3d_comp[1].norm=s3d_comp[1].norm * e2d_norm_optimum(1.e-4);
- Minimum FOM of 5011.2 at norm factor = 0.968792
- - -
hydra> e2d_pars2fom;
...
o Data-Model 10, chi^2 = 4904.13 from 2250 cells; chi^2/cell = 2.17961; sum(D)/sum(M) = 0.990116
o Total FOM = 4904.13
- - -
hydra> e2d_view_resid(10);
The after-fitting residuals, data, and model. FOM = 4904
|  |
|
And take a look at the overall model components in projection for the new cavity parameters:
hydra> s3d_comp[11].rndint=40.0; % adjust for pretty picture. hydra> s3d_proj_comps;Very nice...
- - -
|  |
|
And take a look at the overall model components in projection (left)
and as a 3D solid (ignoring the beta model component):
