HETGS Line Profiles  

ISIS Script for Line Profile Simulations

 

This is the script which generated the figures. It can be run from ISIS via

isis> evalfile("LSF_Script.sl");
or the shortcut,
isis> .load LSF_Script.sl
provided that the script, data, and responses are available in the expected directories. Note that some of the models and fits take several seconds to evaluate. (The Gaussian fits are quite fast and are done first). The script also creates several postscript plots in the working directory. The script also requires the contributed model definitions from aped_fit_models.sl
    %; Time-stamp: <2002-06-06 19:43:25 dph>
    %; MIT Directory: ~dph/d1/CXO_Data/Capella/Sl
    %; File: R_1103.sl
    %; Author: D. Huenemoerder
    %; Original version: 2002.05.08
    %;====================================================================
    % version: 1
    %
    % purpose: Fit a line profile. Simulate many times to examine
    %          statistical deviations from the model profile.
    %          Do for single component line, use gaussian+constant model.
    %          Do for multicomponent line, and use aped plasma model.
    %          In each case, the model is convolved by the rmf*arf.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Reprocessed Capella ObsID 1103; see ../Sh/Proc_1103.sh
    %   de-afterglowed
    %   destreaked
    %   re-did acis_process_events, tg*, ... tgextract,  fullgarf.
    %
    %  NOTE: saturation in core of zero-order, but centroid looks good.
    %  NOTE: osip not so good on S1.  Some Line fluxes 20% low.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    _auto_declare = 1;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%% Load Data, responses, and assign %%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %+
    % Define the data, arfs, rmfs:

    variable f_pha = "./Pha/acisf01103_002N002_pha2.fits";
    variable pha_rows = [1,2,3,4];       % HEG -1,1, then MEG -1,1 rows

    % ARFs in same order as data, for convenience:
    %
    variable f_arfs = "./Arf/acisf01103_002N002" +
           ["HEG_-1", "HEG_1", "MEG_-1", "MEG_1"] +
           "_garf.fits" ;

    % RMFs: HEG, then MEG:
    variable f_rmfs = ["./Ard/acisheg1D1999-07-22rmfN0004.fits",
    		   "./Ard/acismeg1D1999-07-22rmfN0004.fits"];

    % Load the data +-1st orders, HEG and MEG, and associated arfs, and rmfs.
    % Then assign the responses to the data histograms.  Don't save return
    % values.
    %

    () = load_data(f_pha, pha_rows);

    () = load_arf( f_arfs[0] );
    () = load_arf( f_arfs[1] );
    () = load_arf( f_arfs[2] );
    () = load_arf( f_arfs[3] );

    () = load_rmf( f_rmfs[0] );
    () = load_rmf( f_rmfs[1] );

    assign_arf( 1, 1);          % HEG -1
    assign_arf( 2, 2);          % HEG +1
    assign_arf( 3, 3);          % MEG -1
    assign_arf( 4, 4);          % MEG +1

    assign_rmf( 1, [1,2] );      % HEG
    assign_rmf( 2, [3,4] );      % MEG

    flux_corr([1:4]);

    %-


    %%
    % Load the plasma database for use in modeling/line marking.
    % Create a plot window.

    %
    plasma(aped);

    variable p1=open_plot("/xwin"); resize(18, 0.6); erase;

    %%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%% Set up some variables %%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %+

    Fit_Verbose = -1;         % turn off fit output.

    variable n = 24;    %  number of runs to do, 1 w/ real data, 23 fake.

    %%%%%%%%%%%%%%%%%%
    %      n = 3; %%%% FOR QUICK TESTING.....  **** xxxx XXXX DEBUG
    %%%%%%%%%%%%%%%%%

    variable v = Float_Type[n],        % vrad or wavelength array
             f = Float_Type[n],        % flux array
    	 c = Float_Type[n],        % stat array
    	 o = Float_Type[n],        % sigma array
    	 s = Struct_Type;          % fit stat structure

    define Chisq( z )   % utility to reduce chisq from fit stats.
      {
        return z.statistic / (z.num_bins - z.num_variable_params);
      }


    % we will make confidence contours for the observed lines, and overlay
    % contours from one of the simulations.

    variable dcf,dcv, dcs,  % data contour x,y for flux, velocity (or
                            %         wavelength), sigma
             scf, scv, scs, % simulation contour x,y, sigma
    	 dfvmap,  sfvmap, sfsmap;  % data and simulation 2D maps.

    %-


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%  Fit a line parametrically %%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %+
    %
    %%%%%%%  Pick an unblended line  and fit with gaussian + constant
    %

    ignore([1:4]);
    variable xlo=14.165,  xhi=14.235;
    xnotice(2, xlo, xhi);
    xrange(xlo-.02, xhi+.02);
    yrange(0,40);
    label("Wavelength [\\A]", "Flux [counts/bin]", "ObsID 1103/Capella/HEG +1");

    fit_fun( "gauss(1) + poly(1)" );

    set_par( "gauss(1).area", 0.016, 0, 0.0, 0.1 );
    set_par( "gauss(1).center", 14.205, 0, 14.18, 14.23 );
    set_par( "gauss(1).sigma", 0.004, 0, 0.001, 0.008 );

    set_par( "poly(1).a0", 0.01, 0, 0, 0.1 );
    set_par( "poly(1).a1", 0, 1 );
    set_par( "poly(1).a2", 0, 1 );

    () = renorm_counts;

    window(p1);    rplot_counts(2);

    () = fit_counts( &s );
    list_par;

    ()=printf("Conf.sigma: %e, %e\n", conf("gauss(1).sigma"));


    % compute contour map; ad hoc limits.
      dcv = conf_grid("gauss(1).center", 14.201, 14.210, 25);
      dcf = conf_grid("gauss(1).area", 6.e-4, 2.e-3, 25);
      dcs = conf_grid("gauss(1).sigma", 0.0005, 0.01, 25);
      dfvmap = conf_map_counts(dcf, dcv);
      dfsmap = conf_map_counts(dcf, dcs);


    rplot_counts(2);
    freeze( "poly(1).a0" );

    variable p = get_params;    % Copy the parameters into an array of structs.
                                % This will be used to re-instantiate the same
                                % model before the next fakeit.

    variable sp1, sp2, sp3;     % save params vs run;

    %  save the fit norm, vrad, and stat
    %
    f[0] = get_par("gauss(1).area");
    v[0] = get_par("gauss(1).center");
    o[0] = get_par("gauss(1).sigma");
    c[0] = Chisq( s );


    %%% fakeit n-1 times... and display in mosaic w/ the data.
    %%% Save the stat, flux, and centroid from each iteration.

    variable p2=open_plot("/xwin",6,4);resize(8*2.54, 8.5/11. ); erase;

    xrange(xlo-.02, xhi+.02);
    yrange(0,36);
    xnotice(2, xlo, xhi);

    variable pid=dup_plot("1103_14A_HEG_p1_fakeits.ps/cps",p2);
    resize(10*2.54, 8.5/11.);

    window(p2);

    rplot_counts(2);
    window(pid);          % switch to the postscript device...
    rplot_counts(2);      % ... and replot it.
    window(p2);           % go back to screen.

    ignore(2);
    assign_arf(2,5);   % assigning to non-existent histogram index creates one.
    assign_rmf(1,5);
    xnotice(5, xlo, xhi);

    variable i;            % loop over fakeits.....................
    for (i=0;i<n-1;i++) {
      set_params(p);
      fakeit;
      () = fit_counts( &s );
      rplot_counts(5);   window(pid);   rplot_counts(5);   window(p2);
      f[i+1] = get_par("gauss(1).area");
      v[i+1] = get_par("gauss(1).center");
      o[i+1] = get_par("gauss(1).sigma");
      c[i+1] = Chisq( s );
    }


    % save results for summary/comparison of param values and ranges;
    variable f1 = @f;
    variable v1 = @v;
    variable o1 = @o;
    variable sp1 = @p;

    close_plot(pid);
    close_plot(p2);

    % print out results of all fits (could plot these, or
    % do distributions if doing cash statistic and needed distributions
    % for confidence limits)

    vmessage("Moment of fluxes: \n"); print(moment(f));
    vmessage("Moment of center: \n"); print(moment(v));
    vmessage("Moment of sigma: \n"); print(moment(o));
    vmessage("Moment of stat: \n"); print(moment(c));

    % compute contour map for last simulation:
      scv = conf_grid("gauss(1).center", 14.201, 14.21, 25);
      scf = conf_grid("gauss(1).area", 6.e-4, 2.e-3, 25);
      scs = conf_grid("gauss(1).sigma", 0.0005, 0.01, 25);
      sfvmap = conf_map_counts(scf, scv);
      sfsmap = conf_map_counts(scf, scs);

    window(p1);
    label("Flux", "Wavelength", "Conf Maps, data and sim");
    variable line = struct { width, type };
    line.width=2; line.type=2;
    limits;

    plot_conf(dfvmap);  oplot_conf(sfvmap, line);
    connect_points(0); pointstyle(4); oplot(f1, v1, 5 ); connect_points(1);

    pid=dup_plot("1103_14A_HEG_p1_conf_map_1.ps/vcps",p1); resize(6*2.54,1.0);

    plot_conf(dfvmap);  oplot_conf(sfvmap,line);
    connect_points(0); pointstyle(4); oplot(f1, v1, 5 ); connect_points(1);

    label("Flux", "Sigma", "Conf Maps, data and sim");
    plot_conf(dfsmap);  oplot_conf(sfsmap,line);
    connect_points(0); pointstyle(4); oplot(f1, o1, 5 ); connect_points(1);

    close_plot(pid);

    %-

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%% Fit with APED plasma %%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %-
    %  These models take longer to evaluate, so depending on the speed of
    %  the computer, it may be a few seconds per fit, or much longer.
    %
    % define the model functions.
    %
    % aped_fit_models.sl is from dph and defines several plasma functions;
    % the one used here is a single temperature plasma, and has parameters
    % for temperature, radial velocity, turbulent broadening, and
    % normalization.
    %
    % Available from:
    %
    % http://space.mit.edu/CXC/ISIS/contrib/Aped_fit/aped_fit_models.sl
    % (and in the isis contrib directory, with a bug)

    ()=evalfile("aped_fit_models.sl");   % the aped model definitions

    % define a single-temperature plasma model
    %  (list_par will show the params)

    fit_fun("apedx(1)");    % (renamed apedx because of
                            % collision w/ aped variable for database.)

    %%%
    % the following are the result of an ad hoc fit to Ne X 12A:
    %
    set_par("apedx(1).norm", 0.018, 0, 0.0001, 0.08);
    set_par("apedx(1).temp", 10.^6.8, 1, 1.e6, 3.e7);   % a good enough temperature

    % The lines do seem to need more-than-instrumental broadening.  Is
    % thermal broadening  significant for Ne X at 12A in the HEG?

    use_thermal_profile;                          % turn on broadening,
    set_par("apedx(1).vturb", 0, 1, 0.0, 200);    %  but freeze vturb to 0.


    % Since there seems to be a problem with the wavelength scale, we will
    % work around it by letting the radial velocity be free.  (This is OK
    % for single orders, but not for joint fits w/ different offsets.)

    set_par("apedx(1).vrad", 100, 0, 0, 200);   % vrad

    % ignore all but HEG -1 for this test:
    %
    xlo=12.09;   xhi=12.17;
    xnotice(1, xlo, xhi);
    ignore([2:5]);

    yrange(0,80);

    () = renorm_counts;   % adjust the normalization to the data.

    list_data;   % inspect lists
    list_par;

    xrange(xlo-.01, xhi+.01);         % select Ne X + Fe  region
    label("Wavelength [\\A]", "Flux [counts/bin]", "ObsID 1103/Capella/HEG -1");
    plot_data_counts(1,red);
    oplot_model_counts(1,green);

    rplot_counts(1);             % residual plot

    () = fit_counts;

    rplot_counts(1);

    % the profile isn't really as bad as it looks. (but the model isn't
    %   perfect; residuals aren't quite random);
    list_par;


    % compute contour map for data fit
      dcv = conf_grid("apedx(1).vrad", 0.0, 150.0, 20);
      dcf = conf_grid("apedx(1).norm", 0.012, 0.023, 20);
      dfvmap = conf_map_counts(dcf, dcv);


    %  Create some fake data with the same model, and then fit it.
    % do a few times successively.

    assign_arf(1,5);   % assigning to non-existent histogram index creates one.
    assign_rmf(1,5);
    fakeit;                    % uses model to make fake data.
    rplot_counts(5);           % look at it.
    ignore(1);                 % switch fit to histogram 5:
    xnotice(5, xlo, xhi);

    %%%%
    % let's do several fakeits and plot in a m*n window...
    % First go back to histogram 1 and fit, to put that in the window for
    % comparison:

    ignore(5);
    xnotice(1, xlo, xhi);

    %charsize(1.33);  % a little bigger, but not too much to fall out the box.

    p2 = open_plot("/xwin",6,4); resize(8*2.54,8.5/11.); erase;
    xrange(xlo, xhi);
    yrange(0, 80);

    %  set up a ps device.  We will plot to each device.
    %  The ps device will inherit from the xwin if the ranges are set first.
    %
    pid=dup_plot("1103_12A_HEG_m1_fakeits.ps/cps",p2); resize(10.5*2.54, 8.5/11.);
    window(p2);


    () = fit_counts( &s );   % do the fit, and save the statistics structure.
    rplot_counts(1);
    window(pid);
    rplot_counts(1);
    window(p2);

    p = get_params;    % Copy the parameters into an array of structs.
                       % This will be used to re-instantiate the same
                       % model before the next fakeit.

    %  save the fit norm, vrad, and stat:
    %
    f[0] = get_par("apedx(1).norm");
    v[0] = get_par("apedx(1).vrad");
    c[0] = Chisq( s );

    ignore(1);
    xnotice(5, xlo, xhi);

    % Do a series of simulations, each time using the same underlying model.
    %
    for ( i=0; i<n-1; i++ ) {
      set_params( p );
      fakeit;
      () = fit_counts( &s );
      rplot_counts(5);  window(pid); rplot_counts(5); window(p2);
      f[i+1] = get_par("apedx(1).norm");
      v[i+1] = get_par("apedx(1).vrad");
      c[i+1] = Chisq( s );
    }
    close_plot(pid);

    f2 = @f;
    v2 = @v;
    sp2 = @p;

    vmessage("Moment of fluxes: \n"); print(moment(f));
    vmessage("Moment of vrad: \n"); print(moment(v));
    vmessage("Moment of stat: \n"); print(moment(c));

    % compute contour map for last simulation:
      scv = conf_grid("apedx(1).vrad", 0.0, 150.0, 20);
      scf = conf_grid("apedx(1).norm", 0.012, 0.023, 20);
      sfvmap = conf_map_counts(scf, scv);

    window(p1);
    label("Flux", "Radial Velocity", "Conf Maps, data and sim");
    line.width=2; line.type=2;
    limits;

    plot_conf(dfvmap);  oplot_conf(sfvmap, line);
    connect_points(0); pointstyle(4); oplot(f2, v2, 5 );connect_points(1);

    pid=dup_plot("1103_12A_HEG_m1_conf_map_2.ps/vcps",p1); resize(6*2.54,1.0);
    plot_conf(dfvmap);  oplot_conf(sfvmap,line);
    connect_points(0); pointstyle(4); oplot(f2, v2, 5 );connect_points(1);
    close_plot(pid);

    %-

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% plot the flux model, to see lines:

    window(p1);
    ignore(5);
    xnotice(1, xlo, xhi);
    xrange(xlo, xhi);
    yrange(0, 6.e-4);
    label("Wavelength [\\A]", "Flux [counts/bin]", "ObsID 1103/Capella/HEG -1");

    pid = dup_plot("1103_12A_model_flux_H_m1.ps/cps", p1);  % set up ps
    window(p1);                                   % ...but select xwin

    set_par("apedx(1).vrad", 0 );  % turn off shift

    use_delta_profile;
    eval_counts;
    plot_model_flux(1,2);
    plot_group(brightest(10,where(wl(12.1,12.17))),1);

    window(pid);
    plot_model_flux(1,2);
    plot_group(brightest(10,where(wl(12.1,12.17))),1);
    window(p1);

    use_thermal_profile;
    eval_counts;
    oplot_model_flux(1,3);

    window(pid);
    oplot_model_flux(1,3);

    close_plot(pid); window(p1);
    close_plot(p2);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%% Repeat, but without thermal broadening, to see resids:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %-

    use_delta_profile;
    set_params( p );

    xlo=12.09;   xhi=12.17;
    xnotice(1, xlo, xhi);
    ignore([2:5]);
    yrange(0,80);

    p2 = open_plot("/xwin",6,4); resize(8*2.54,8.5/11.); erase;
    xrange(xlo, xhi);
    yrange(0, 80);

    pid=dup_plot("1103_12A_HEG_m1_fakeits_delta.ps/cps",p2); resize(10.5*2.54, 8.5/11.);
    window(p2);

    () = fit_counts( &s );   % do the fit, and save the statistics structure.
    rplot_counts(1);
    window(pid);
    rplot_counts(1);
    window(p2);

    p = get_params;    % Copy the parameters into an array of structs.
                       % This will be used to re-instantiate the same
                       % model before the next fakeit.

    %  save the fit norm, vrad, and stat:
    %
    f[0] = get_par("apedx(1).norm");
    v[0] = get_par("apedx(1).vrad");
    c[0] = Chisq( s );

    % compute contour map for data fit
      dcv = conf_grid("apedx(1).vrad", 0.0, 150.0, 20);
      dcf = conf_grid("apedx(1).norm", 0.012, 0.023, 20);
      dfvmap = conf_map_counts(dcf, dcv);

    ignore(1);
    xnotice(5, xlo, xhi);

    % Do a series of simulations, each time using the same underlying model.
    %
    for ( i=0; i<n-1; i++ ) {
      set_params( p );
      fakeit;
      () = fit_counts( &s );
      rplot_counts(5);  window(pid); rplot_counts(5); window(p2);
      f[i+1] = get_par("apedx(1).norm");
      v[i+1] = get_par("apedx(1).vrad");
      c[i+1] = Chisq( s );
    }
    close_plot(pid);

    v3 = @v;
    f3 = @f;
    sp3 = @p;

    vmessage("Moment of fluxes: \n"); print(moment(f));
    vmessage("Moment of dlambda: \n"); print(moment(v));
    vmessage("Moment of stat: \n"); print(moment(c));

    % compute contour map for last simulation:
      scv = conf_grid("apedx(1).vrad", 0.0, 150.0, 20);
      scf = conf_grid("apedx(1).norm", 0.012, 0.023, 15);
      sfvmap = conf_map_counts(dcf, dcv);

    window(p1);
    label("Flux", "Radial Velocity", "Conf Maps, data and sim");
    line.width=2; line.type=2;
    limits;

    plot_conf(dfvmap);  oplot_conf(sfvmap, line);
    connect_points(0); pointstyle(4); oplot(f3, v3, 5 ); connect_points(1);
    pid=dup_plot("1103_12A_HEG_m1_conf_map_3.ps/vcps",p1); resize(6*2.54,1.0);
    plot_conf(dfvmap);  oplot_conf(sfvmap,line);
    connect_points(0); pointstyle(4); oplot(f3, v3, 5 ); connect_points(1);
    close_plot(pid);

    %-

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Fit_Verbose = 0;


Previous: Figures Next: Data


This page was last updated Jun 6, 2002 by David P. Huenemoerder. To comment on it or the material presented here, send email to dph@space.mit.edu.
Valid HTML 4.01! Made with JED Viewable With Any Browser