## 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.

- The data are read in twice, once for X-Y coord.s and once for RADIUS-ENERGY coord.s.

- Two obsids are being combined to make each single data set - see the code/comments in

- The X,Y offsets on input for the obsids were adjusted to put the AGN at ~ 0,0 in the e2d coord.s. We'll see that the best-fit cluster emission is somewhat offset from this.

- The "zo" means zeroth-order and the data are put in id's 10 and 11. This is a hold-over from the practice with grating datasets to put the dispersed data in datasets 1,2,3,4, and the zeroth-order imaging data in dataset 10 (and 11).

hydra>Look at the datasets....source ha_zodataReading 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 A

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>The model-defining file,.source ha_model_simpleSpectrum 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 ...

--- 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} - - -

NOTE: The following was done with an incorrect definition of beta in the model routine.

The model has been updated and now beta values of 0.433 and 0.411 correspond to the previous values of 0.6 and 0.578.

Keep this in mind if duplicating the following.

Evaluate the source model... (This was already done in the ha_model_simple.sl file - this would be used if some parameters were manually changed.)

hydra>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 thes3d_update;o check/update foreground NH ... o check/update spectra ... o check/update the component geometry ...

hydra>Show the spectra and NH curve...s3d_view_comps;hydra>s3d_proj_comps;

hydra>s3d_plot_spectrum(19);% 19 is the AGN spectrum Graphics device/type (? to see list, default /XSERVE):CRhydra>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>Look at the Data, Model and Residual images...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 - - -

hydra>At Left: Color-coded residual image - see print-out above for color coding info. At Right: The observed data, +/-150" range.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 ]

Below-left is the simulated data from a +/-100" cube.

- - -Re-adjust the viewing (and FOM calculating) ranges to extend to only 100" to be within the simulation size.

hydra>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 !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]);

hydra>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.)e2d_plot_ks(10,1);hydra>e2d_plot_ks(10,2);hydra>e2d_plot_ks(11,1);hydra>e2d_plot_ks(11,2);

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>and calculate this adjustment more accurately with: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

hydra>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:norm_adjust = e2d_norm_optimum(1.e-4);- Minimum FOM of 40316.4 at norm factor = 1.20653 - - -

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 = 39693OK, that's better !- - -

The present pretty-good values were set by manual adjustment... Will fitting find much different ? We'll see... But in any case we'll see how the fitting works :-)

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>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.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 ====================

- - -Now, set e2d_update_pars to scan the beta value (after putting the betax value back to manual value):

hydra>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.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 ====================

- - -

- Set the dimension of the variable for the fit parameters and their scale values, aka "smand" ("should make a noticable difference") values.

- Set the starting location in parameter space.

- Set the smand values.

- And, finally, define the e2d_update_pars that is appropriate for cg fitting for these parameters. As in the simple scan, it defines a 1D scan but now in an arbitrary direction, cg_dir, in the parameter space. The smand values serve as the unit-distances for each direction, performing a simple "pre-conditioning" to reduce the ellipticity of the FOM contours in the search space.

hydra>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...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]; }

hydra>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: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.

**
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

beta-center (x,y) = (-1.44", 1.28") and beta = 0.573

with FOM ~ 38000.
Note, however, that the gradients (dFOM/dParameter) are very
different for the beta parameter than for the beta-x,y parameters,
as shown after the first gradient evaluation:

-1*Gradient(p0) = -2858.02, 2734.2, 45107.6

Better match these by adjusting the smand values,
especially by reducing the size of the beta smand value quite a bit,
and continue the fit starting at the place it left off:

hydra>OK, now after ~ 80 minutes we get: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 - switch

In this case the fitting stopped after evaluating the local derivatives at the 5th iteration because all of the dFOM/dParameter values were within the noise level of 0 - hence the routine has no preferred direction to search and it stops. Note that it almost stopped after the 2nd iteration except that the FOM had changed a bit more than the stopping criteria. For this fitting, the final values are:

beta-center (x,y) = (-2.40", 1.98") and beta = 0.579

with FOM ~ 37400.

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:

beta-center (x,y) = (-2.47", 2.31") and beta = 0.578

also with FOM ~ 37400. The data-model residuals
for this final case can be viewed by doing:

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>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..source ha_model_cavityhydra>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;

- - -

hydra>One cavity in the southwest: starting residuals, data, and model..source ha_zodatahydra>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 ]

At Left: Color-coded residual image. At Right: The observed data.

Below-left is the simulated data which includes a cavity.

- - -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...

Ouch! - that was too large a scale - it sent the radius negative... Try again with more modest smand's values of [10.0,10.0,3.0] :

hydra>Here's the log of what the fitting did for about 3 hours: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);

Looks like it did a good job of consistantly moving to lower FOMs and the parameter values are converging after the 8 iterations to something like: (x,y) about (+40., -54.) and a radius of 33". Let's fix the parameters to these values and redo the image and residuals we did above:

hydra>The after-fitting residuals, data, and model.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);

At Left: Color-coded residual image. At Right: The observed data.

Below-left is the simulated data which includes the fit cavity.

And take a look at the overall model components in projection for the new cavity parameters:

hydra>Very nice...s3d_comp[11].rndint=40.0;% adjust for pretty picture. hydra>s3d_proj_comps;

- - -

Now fit other cavities as well using the commands and recording progress in the file

to get residuals, data, model images:

And take a look at the overall model components in projection (left)
and as a 3D solid (ignoring the beta model component):