Hydra A - Event-2D Demo



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


Reading-in and Looking-at the Data

Load in the user-defined data sets. Some things to note:
- 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 ha_zodata.sl for details.
- 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> .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 A

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

Reading-in and Looking-at a Model

Load in and list the user-defined model which consists of Spectra and geometric Components...
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}
 - - -


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

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

 - - -

Comparing Data and Model

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

 - - -

Adjusting the Model Norm

OK, let's concentrate on the X,Y dataset, 10 by ignoring dataset 11 and checking the FOM. Then adjust the norm by scanning through some scale factors for it and seeing the FOM vs factor plot:
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 !
 - - -

Almost Fitting: Scanning Parameters

Now let's do some fitting for the best X,Y location and beta parameter for the simple model!
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> 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.

 - - -

Fitting with Congugate Gradient routine

The routine e2d_cg_fit does a nonlinear conjugate gradient fit using the e2d_scan_par routine to search/fit along a given 1D direction. To set things up to fit the three parameters - betax, betay, and beta - there are just a few things to setup; it's not too tricky:
- 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> 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...
"And the monkey flips the switch !"
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
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> 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 
OK, now after ~ 80 minutes we get:
cg_log_Mon_Nov_12_16_54_38_2007.html
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... :-)

Fitting a Cavity

Use a model that includes spherical cavities in the beta emission and, for demonstration, define and have the fitting adjust the parameters of a single cavity... This model also includes a wisp of emission in the NE, seen in the image below.
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
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> 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:
cg_log_Tue_Nov_13_10_56_33_2007.html
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> 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
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> s3d_comp[11].rndint=40.0;   % adjust for pretty picture.
hydra> s3d_proj_comps;

- - -
Very nice...
Now fit other cavities as well using the commands and recording progress in the file analysis_cavityJ.sl...
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):