with Kevin, 5 August 2005 - 

Going through procedure to fit a line in the spectrum with different models in xspec, generate a steppar log, use Kevin's programs to plot the data plus model, and determine confidence limits. 

*Our test problem is the line (complex) near 16.00 Angstroms in zeta Ori

First, we identify the line(s) and look for blends: 

http://cxc.harvard.edu/atomdb/WebGUIDE/index.html

There are two O VIII transitions (Lyman beta):  16.0055, 16.0067 (2:1 intensity ratio, so the mean wavelength is 16.0059; note that the splitting between the two lines is negligible - 22 km/sec)

Blends? 

Fe XVIII at 16.0040 (formation temperature is twice that of the O VIII) NOTE: wavelength is nearly identical, so we can ignore this separate component

Fe XVIII at 16.0710 (also a formation temp twice that of O VIII; relative intensity is 0.84) - so we should model this. 

note that there's an Fe XIX line at 16.1100 (formation temperature is about three times higher than O VIII, but rel. intensity is 0.46)

*we will fit a two component model to this feature: one component (O VIII doublet plus Fe XVIII) at 16.0059 and the other (Fe XVIII) at 16.0710

Now, we will look at the spectrum in xspec: 

working directory with zeta Ori data, etc.: 

/home/cohen/xspec/2005/

[cohen@ruby 2005]$ xspec

               Xspec 11.3.1    14:38:50 05-Aug-2005

 For documentation, notes, and fixes see  http://xspec.gsfc.nasa.gov/

 Plot device not set, use "cpd" to set it

 Type "help" or "?" for further information
XSPEC>

XSPEC>@loaddata  		!Maurice's script for reading in the data, responses, etc. (loaddata.xcm)

Aside: Kevin also set up a small script: 

[cohen@ruby 2005]$ more setup.xcm
setplot wave
query yes
statistic cstat
cpd /xw

but he says that the cpd command must be typed in by hand

XSPEC>@setup

!XSPEC>  setplot wave;

!XSPEC>  query yes;
Querying disabled - assuming answer is yes

!XSPEC>  statistic cstat;
Warning : this statistic is only valid for Poisson data

XSPEC>cpd /xw &			!open a xwindow

Note: to look at and analyze only one order, edit the loaddata script

XSPEC>plot				!look at all the data

Now, choosing the wavelength range to look at (and fit): 

XSPEC>ignore 0.-15.8
 ignoring channels  5233 -  8192 in dataset     1
 ignoring channels  5233 -  8192 in dataset     2
Note that RESPFILE keyword in ARF is grid(acisf00610MEG-1_grmf.fits[cols ENERG_LO, ENERG_HI])
Note that RESPFILE keyword in ARF is grid(acisf00610MEG1_grmf.fits[cols ENERG_LO, ENERG_HI])
XSPEC>ignore 16.2-** 
 ignoring channels     1 -  5153 in dataset     1
 ignoring channels     1 -  5153 in dataset     2
Note that RESPFILE keyword in ARF is grid(acisf00610MEG-1_grmf.fits[cols ENERG_LO, ENERG_HI])
Note that RESPFILE keyword in ARF is grid(acisf00610MEG1_grmf.fits[cols ENERG_LO, ENERG_HI])

So, this has 15.8 to 16.2 as the data range

(Aside: the x-axis label bug is preventing zoomed in views from showing labelled x-axes)

XSPEC>plot

Looking at the data, I see no obvious evidence of the blends near 16.07 and 16.11, so first, I'll fit just a single component model: 

1) Fixed wavelength Gaussian (plus powerlaw): 

XSPEC>model lgauss + power
  Model:  lgauss<1> + powerlaw<2>
Input parameter value, delta, min, bot, top, and max values for ...
               18.97      -0.1         1         1       200       200
1:lgauss:waveleng>16.0059
                0.01      -0.1         0         0         1         1
2:lgauss:sigma>5e-2
                   1      0.01         0         0     1E+24     1E+24
3:lgauss:norm>5.e-4
                   1      0.01        -3        -2         9        10
4:powerlaw:PhoIndex>2.
                   1      0.01         0         0     1E+24     1E+24
5:powerlaw:norm>5.e-4
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  lgauss<1> + powerlaw<2>
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   lgauss     waveleng A          16.0059     frozen
    2    2    1   lgauss     sigma    A         5.000000E-02 frozen
    3    3    1   lgauss     norm               5.000000E-04 +/-   0.00000
    4    4    2   powerlaw   PhoIndex            2.00000     +/-   0.00000
    5    5    2   powerlaw   norm               5.000000E-04 +/-   0.00000
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 C-statistic =      490.6682     using   158 PHA bins.

XSPEC>thaw 2
 Number of variable fit parameters =    4
XSPEC>freeze 4
 Number of variable fit parameters =    3

XSPEC>show 
...
  Model:  lgauss<1> + powerlaw<2>
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  lgauss<1> + powerlaw<2>
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   lgauss     waveleng A          16.0059     frozen
    2    2    1   lgauss     sigma    A         5.000000E-02 +/-   0.00000
    3    3    1   lgauss     norm               5.000000E-04 +/-   0.00000
    4    4    2   powerlaw   PhoIndex            2.00000     frozen
    5    5    2   powerlaw   norm               5.000000E-04 +/-   0.00000
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 C-statistic =      490.6682     using   158 PHA bins.

XSPEC>plot

after typing "fit" - 

...
   178.330      0     3.5746E-02  1.6125E-04  3.1979E-03
 ---------------------------------------------------------------------------
 Variances and Principal axes :
               2     3     5
 1.85E-10 |  0.00  1.00  0.00
 4.26E-03 | -1.00  0.00  0.03
 2.91E-07 | -0.03  0.00 -1.00
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  lgauss<1> + powerlaw<2>
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   lgauss     waveleng A          16.0059     frozen
    2    2    1   lgauss     sigma    A         3.574598E-02 +/-  0.652507E-01
    3    3    1   lgauss     norm               1.612459E-04 +/-  0.101153E-03
    4    4    2   powerlaw   PhoIndex            2.00000     frozen
    5    5    2   powerlaw   norm               3.197896E-03 +/-  0.184450E-02
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 C-statistic =      178.3300     using   158 PHA bins.

XSPEC>plot

XSPEC>goodness 100
 61.00% of realizations have a fit statistic <      178.3

This fit is formally good.  But, I'm wondering whether there's a better set of parameters.  Type "fit" again - it improved a little: 

XSPEC>fit
 C-statistic    Lvl  Fit param # 2     3           5
   178.324      0     3.5733E-02  1.6109E-04  3.2003E-03
 ---------------------------------------------------------------------------
 Variances and Principal axes :
               2     3     5
 1.84E-10 |  0.00  1.00  0.00
 4.27E-03 | -1.00  0.00  0.03
 2.91E-07 | -0.03  0.00 -1.00
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  lgauss<1> + powerlaw<2>
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   lgauss     waveleng A          16.0059     frozen
    2    2    1   lgauss     sigma    A         3.573310E-02 +/-  0.653112E-01
    3    3    1   lgauss     norm               1.610895E-04 +/-  0.101147E-03
    4    4    2   powerlaw   PhoIndex            2.00000     frozen
    5    5    2   powerlaw   norm               3.200253E-03 +/-  0.184456E-02
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 C-statistic =      178.3238     using   158 PHA bins.

But, we can also use steppar to find a better fit: 

XSPEC>steppar 2 2.e-2 5.e-2 
 C-statistic        Delta   Mod param
              C-statistic         2
     215.84        37.520  2.000000E-02
     204.27        25.945  2.300000E-02
     194.29        15.961  2.600000E-02
...

Whoops, I forgot the last parameter, which is the number of steps; it appears to be using some default number, maybe 10

...
     186.48         8.157  2.900000E-02
     181.27         2.945  3.200000E-02
     178.40         0.081  3.500000E-02
     177.78        -0.539  3.800000E-02
     178.72         0.393  4.100000E-02
     181.14         2.820  4.400000E-02
     184.49         6.169  4.700000E-02
     188.52        10.196  5.000000E-02
XSPEC>

See, it found a better fit for sigma=3.8e-2

XSPEC>newpar 2 3.8e-2
    3 variable fit parameters
 C-statistic =      177.7794     using   158 PHA bins.

*Kevin says, sometimes after changing a parameter value with the newpar command, typing fit won't budge it.  So you freeze it, then you fit, then you thaw, then re-fit. 

XSPEC>goodness 1000
 56.90% of realizations have a fit statistic <      177.8

OK, this fit is pretty good.  But (in a minute, we'll fit thawed G. and windprofile) - first, though - let's fit two Gaussians with fixed wavelengths.  

XSPEC>model lgauss + lgauss + pow
  Model:  lgauss<1> + lgauss<2> + powerlaw<3>
Input parameter value, delta, min, bot, top, and max values for ...
               18.97      -0.1         1         1       200       200
1:lgauss:waveleng>16.0059
                0.01      -0.1         0         0         1         1
2:lgauss:sigma>4.e-2
                   1      0.01         0         0     1E+24     1E+24
3:lgauss:norm>1.6e-4
               18.97      -0.1         1         1       200       200
4:lgauss:waveleng>16.0710
                0.01      -0.1         0         0         1         1
5:lgauss:sigma>4.e-2
                   1      0.01         0         0     1E+24     1E+24
6:lgauss:norm>1.e-4
                   1      0.01        -3        -2         9        10
7:powerlaw:PhoIndex>2.
                   1      0.01         0         0     1E+24     1E+24
8:powerlaw:norm>3.e-3

XSPEC>thaw 2 
 Number of variable fit parameters =    5
XSPEC>thaw 5
 Number of variable fit parameters =    6
XSPEC>freeze 7
 Number of variable fit parameters =    5
XSPEC>show
...

Model:  lgauss<1> + lgauss<2> + powerlaw<3>
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  lgauss<1> + lgauss<2> + powerlaw<3>
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   lgauss     waveleng A          16.0059     frozen
    2    2    1   lgauss     sigma    A         4.000000E-02 +/-   0.00000
    3    3    1   lgauss     norm               1.600000E-04 +/-   0.00000
    4    4    2   lgauss     waveleng A          16.0710     frozen
    5    5    2   lgauss     sigma    A         4.000000E-02 +/-   0.00000
    6    6    2   lgauss     norm               1.000000E-04 +/-   0.00000
    7    7    3   powerlaw   PhoIndex            2.00000     frozen
    8    8    3   powerlaw   norm               3.000000E-03 +/-   0.00000
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 C-statistic =      271.5660     using   158 PHA bins.

XSPEC>newpar 5 = 2
Equating parameter :sigma to parameter lgauss:sigma  * 1
    4 variable fit parameters
 C-statistic =      271.5660     using   158 PHA bins.

This last command tied the two model sigmas together

  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  lgauss<1> + lgauss<2> + powerlaw<3>
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   lgauss     waveleng A          16.0059     frozen
    2    2    1   lgauss     sigma    A         3.128324E-02 +/-  0.721712E-01
    3    3    1   lgauss     norm               1.307371E-04 +/-  0.905009E-04
    4    4    2   lgauss     waveleng A          16.0710     frozen
    5    2    2   lgauss     sigma    A         3.128324E-02 = par   2
    6    5    2   lgauss     norm               1.426705E-05 +/-  0.112077E-04
    7    6    3   powerlaw   PhoIndex            2.00000     frozen
    8    7    3   powerlaw   norm               3.580004E-03 +/-  0.180438E-02
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 C-statistic =      186.4230     using   158 PHA bins.

Note: tied sigmas and 2nd norm only about 1/10 of first; but also note,  this fit is worse than the one G fit. 

We'll steppar to find a better fit: 

XSPEC>steppar 2 2.5e-2 4.e-2 8
 C-statistic        Delta   Mod param
              C-statistic         2
     196.12         9.696  2.500000E-02
     191.28         4.852  2.687500E-02
     187.00         0.575  2.875000E-02
     183.65        -2.776  3.062500E-02
     180.55        -5.870  3.250000E-02
     178.79        -7.635  3.437500E-02
     177.92        -8.504  3.625000E-02
     177.86        -8.567  3.812500E-02
     178.49        -7.936  4.000000E-02

So, there are much better fits with bigger sigmas (but I bet, smaller norm2s)

XSPEC>newpar 2 3.81e-2

XSPEC>freeze 2
 Number of variable fit parameters =    3
XSPEC>fit
 C-statistic    Lvl  Fit param # 3     5           7
   182.870     -3     1.6263E-04  1.8708E-06  3.6550E-03
   181.410     -4     1.3355E-04  8.4290E-07  3.5302E-03
   181.308     -1     1.5846E-04  2.5050E-07  3.5807E-03
   180.628     -2     1.3889E-04  7.2923E-08  3.5628E-03
   180.293      0     1.4675E-04  1.5129E-08  3.5643E-03
   180.271      0     1.4766E-04  6.3033E-09  3.5353E-03
   180.263      0     1.4789E-04  1.9689E-09  3.5230E-03
 ---------------------------------------------------------------------------
 Variances and Principal axes :
               3     5     7
 2.18E-10 | -0.29 -0.96 -0.01
 7.67E-10 |  0.96 -0.29  0.01
 4.07E-07 |  0.01  0.01 -1.00
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  lgauss<1> + lgauss<2> + powerlaw<3>
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   lgauss     waveleng A          16.0059     frozen
    2    2    1   lgauss     sigma    A         3.275010E-02 frozen
    3    3    1   lgauss     norm               1.478855E-04 +/-  0.278857E-04
    4    4    2   lgauss     waveleng A          16.0710     frozen
    5    2    2   lgauss     sigma    A         3.275010E-02 = par   2
    6    5    2   lgauss     norm               1.968890E-09 +/-  0.174657E-04
    7    6    3   powerlaw   PhoIndex            2.00000     frozen
    8    7    3   powerlaw   norm               3.523039E-03 +/-  0.637868E-03
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 C-statistic =      180.2632     using   158 PHA bins.

... 
the C stat got a little lower, but only by sending the normalization of the second component close to zero.

*Just for the sake of seeing how to do it, let's do a 2 Gaussian thawed but tied wavelength 

XSPEC>newpar 4 = 1+0.0651
Equating parameter :waveleng to parameter lgauss:waveleng  * 1 +0.0651
    5 variable fit parameters
 C-statistic =      179.0650     using   158 PHA bins.

*It found a better fit: 

  ---------------------------------------------------------------------------
  Model:  lgauss<1> + lgauss<2> + powerlaw<3>
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   lgauss     waveleng A          15.9858     +/-  0.624502E-02
    2    2    1   lgauss     sigma    A         2.811126E-02 +/-  0.682481E-01
    3    3    1   lgauss     norm               1.259576E-04 +/-  0.819945E-04
    4    1    2   lgauss     waveleng A          16.0509     = par   1 +0.651000E-01
    5    2    2   lgauss     sigma    A         2.811126E-02 = par   2
    6    4    2   lgauss     norm               3.302195E-05 +/-  0.209574E-04
    7    5    3   powerlaw   PhoIndex            2.00000     frozen
    8    6    3   powerlaw   norm               3.120462E-03 +/-  0.187080E-02
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 C-statistic =      173.5380     using   158 PHA bins.

XSPEC>goodness 100
 44.00% of realizations have a fit statistic <      173.5

Scientifically, though, the conclusion is no real need for the second model component. 

*Now, fit a single Gaussian with a thawed wavelength: 

  ---------------------------------------------------------------------------
  Model:  lgauss<1> + powerlaw<2>
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   lgauss     waveleng A          16.0007     +/-  0.640549E-02
    2    2    1   lgauss     sigma    A         3.569465E-02 +/-  0.738643E-01
    3    3    1   lgauss     norm               1.507715E-04 +/-  0.106895E-03
    4    4    2   powerlaw   PhoIndex            2.00000     frozen
    5    5    2   powerlaw   norm               3.396277E-03 +/-  0.196290E-02
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 C-statistic =      174.6505     using   158 PHA bins.

Can it do better? 

XSPEC>steppar 1 15.99 16.01 12
 C-statistic        Delta   Mod param
              C-statistic         1
     180.94         6.288   15.9900
     178.80         4.145   15.9917
     177.08         2.427   15.9933
     175.79         1.136   15.9950
     174.93         0.278   15.9967
     174.50        -0.146   15.9983
     174.52        -0.135   16.0000
     174.96         0.313   16.0017
     175.85         1.196   16.0033
     177.16         2.512   16.0050
     178.91         4.259   16.0067
     181.09         6.435   16.0083
     183.61         8.961   16.0100

after much stepparing - 

 ---------------------------------------------------------------------------
  Model:  lgauss<1> + powerlaw<2>
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   lgauss     waveleng A          15.9991     +/-  0.631684E-02
    2    2    1   lgauss     sigma    A         3.565472E-02 +/-  0.739812E-01
    3    3    1   lgauss     norm               1.511182E-04 +/-  0.107403E-03
    4    4    2   powerlaw   PhoIndex            2.00000     frozen
    5    5    2   powerlaw   norm               3.381695E-03 +/-  0.197691E-02
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 C-statistic =      174.4390     using   158 PHA bins.

XSPEC>goodness 1000
 47.50% of realizations have a fit statistic <      174.4

Next, fitting a windprofile: 

XSPEC>model windprof + pow
  Model:  windprof<1> + powerlaw<2>
Input parameter value, delta, min, bot, top, and max values for ...
                   0      -0.1     -0.99     -0.99         5         5
1:windprof:q>0
                   1      -0.1         0         0       100       100
2:windprof:taustar>0.3
                 0.5      -0.1         0         0      0.99      0.99
3:windprof:umax>
                2485      -0.1       500       500     5E+03     5E+03
4:windprof:velocity>1860
                   1      -0.1         1         1         1         1
5:windprof:beta>
               24.78      -0.1         1         1       200       200
6:windprof:waveleng>16.0059
                   0      -0.1         0         0         1         1
7:windprof:verbose>
                   1      0.01         0         0     1E+24     1E+24
8:windprof:norm>2.e-4
                   1      0.01        -3        -2         9        10
9:powerlaw:PhoIndex>2.
                   1      0.01         0         0     1E+24     1E+24
10:powerlaw:norm>3.e-3
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  windprof<1> + powerlaw<2>
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   windprof   q                   0.00000     frozen
    2    2    1   windprof   taustar            0.300000     frozen
    3    3    1   windprof   umax               0.500000     frozen
    4    4    1   windprof   velocity km/s       1860.00     frozen
    5    5    1   windprof   beta                1.00000     frozen
    6    6    1   windprof   waveleng A          16.0059     frozen
    7    7    1   windprof   verbose             0.00000     frozen
    8    8    1   windprof   norm               2.000000E-04 +/-   0.00000
    9    9    2   powerlaw   PhoIndex            2.00000     +/-   0.00000
   10   10    2   powerlaw   norm               3.000000E-03 +/-   0.00000
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 C-statistic =      185.3912     using   158 PHA bins.
XSPEC>plot

after thawing and freezing appropriate params: 

  ---------------------------------------------------------------------------
  Model:  windprof<1> + powerlaw<2>
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   windprof   q                 -0.228432     +/-  0.409098
    2    2    1   windprof   taustar            0.255593     +/-  0.226670
    3    3    1   windprof   umax               0.643415     +/-  0.796848E-01
    4    4    1   windprof   velocity km/s       1860.00     frozen
    5    5    1   windprof   beta                1.00000     frozen
    6    6    1   windprof   waveleng A          16.0059     frozen
    7    7    1   windprof   verbose             0.00000     frozen
    8    8    1   windprof   norm               1.623932E-04 +/-  0.277625E-04
    9    9    2   powerlaw   PhoIndex            2.00000     frozen
   10   10    2   powerlaw   norm               3.233301E-03 +/-  0.619411E-03
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 C-statistic =      174.7456     using   158 PHA bins.

  ---------------------------------------------------------------------------
  Model:  windprof<1> + powerlaw<2>
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   windprof   q                 -0.302494     +/-  0.343511
    2    2    1   windprof   taustar            0.270245     +/-  0.223252
    3    3    1   windprof   umax               0.662971     +/-  0.746429E-01
    4    4    1   windprof   velocity km/s       1860.00     frozen
    5    5    1   windprof   beta                1.00000     frozen
    6    6    1   windprof   waveleng A          16.0059     frozen
    7    7    1   windprof   verbose             0.00000     frozen
    8    8    1   windprof   norm               1.550820E-04 +/-  0.286634E-04
    9    9    2   powerlaw   PhoIndex            2.00000     frozen
   10   10    2   powerlaw   norm               3.242555E-03 +/-  0.619393E-03
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 C-statistic =      174.2874     using   158 PHA bins.

After much stepparing we have our best-fit model

XSPEC>iplot 
PLT> wdata windprof_160059 
PLT> exit

XSPEC>ls
acisf00610MEG_-1_garf.fits    acisf00610N002_pha2.fits  notes_may05.txt   test.txt
acisf00610MEG_1_garf.fits     cycle7                    notes_may05.txt~  windprof_160059.qdp
acisf00610MEG-1_grmf.fits     david_test2.txt           setup.xcm         xautosav.xcm
acisf00610MEG1_grmf.fits      idl                       testii.dat
acisf00610N002_bkg_pha2.fits  loaddata.xcm              testiii.dat

So, the windprof_160059.qdp file has the data and model in it: 

[cohen@ruby 2005]$ more windprof_160059.qdp 
READ SERR 1 2
@windprof_160059.pco
!
16.1926575 0.00250021904 0.00335365627 0.00780129246 0.0034525746
16.1876583 0.00249993498 0 0.00627261493 0.00345565798
16.1826572 0.00250027911 0.00670715142 0.00892460905 0.00345805474
16.1776581 0.00249999249 -0.000834561419 0.00628633518 0.00346148829
16.172657 0.00249970471 0.00335434615 0.00780289713 0.00346511835
16.1676579 0.00250004441 0 0.00627234019 0.00346806715
16.1626568 0.00250038202 0.00335343787 0.00780078396 0.00347111165
16.1576576 0.0024994628 0.00293730199 0.0078036529 0.0034759345
16.1526585 0.00250042463 0 0.00627138652 0.00347839156
16.1476574 0.00250013033 0 0.00627212506 0.00348342722
16.1426563 0.0024998351 0.00335417152 0.00780249108 0.00348989316
16.1376572 0.00250016456 0.00335372938 0.00780146243 0.00349717611
16.132658 0.00249986676 0 0.0062727863 0.00350713171
16.1276569 0.00250019296 0.00335369119 0.00780137349 0.00351925753
16.1226578 0.00249989261 0 0.0062727211 0.00354200532
16.1176567 0.00250021555 0.00670732185 0.00892483629 0.00359077868
16.1126575 0.00249991263 0.00335406745 0.00780224847 0.00370708387

columns: 

1: wavelength
2: bin half-width
3: data (counts/sec/angstrom)
4: errors (same units)
5: model (same units)

*Remaining tasks: 

a) plot the qdp file contents using Kevin's IDL program
b) run a steppar log so that we can
c) use Kevin's other program for calculating confidence limits and making contour plots

---
b) To run steppar, first estimate the range of each of the three parameters, using a single, coarse steppar run: 

XSPEC>steppar 1 -0.9 1.2 14
 C-statistic        Delta   Mod param
              C-statistic         1
     187.15        12.865 -0.900000
     178.71         4.422 -0.750000
     175.63         1.345 -0.600000
     174.55         0.261 -0.450000
     174.29         0.000 -0.300000
     174.39         0.098 -0.150000
     174.57         0.280  3.405980E-08
     175.53         1.239  0.150000
     174.87         0.580  0.300000
     175.01         0.724  0.450000
     175.04         0.757  0.600000
     176.33         2.046  0.750000
     176.59         2.300  0.900000
     176.80         2.513   1.05000
     178.04         3.757   1.20000

So, for q, maybe we choose a range of -0.75 to 1.25

XSPEC>steppar 2 0.0 1.2  7
 C-statistic        Delta   Mod param
              C-statistic         2
     177.33         3.044   0.00000
     174.55         0.261  0.171429
     174.39         0.099  0.342857
     175.19         0.905  0.514286
     176.47         2.179  0.685714
     178.00         3.714  0.857143
     179.63         5.341   1.02857
     181.35         7.064   1.20000

And for taustar, we'll use 0 to 1.0

XSPEC>steppar 3 0.3 0.99 8
 C-statistic        Delta   Mod param
              C-statistic         3
     179.41         5.120  0.300000
     176.34         2.055  0.386250
     175.19         0.904  0.472500
     174.73         0.446  0.558750
     174.35         0.058  0.645000
     174.88         0.591  0.731250
     177.76         3.474  0.817500
     182.57         8.278  0.903750
     184.77        10.479  0.990000

for umax, use 0.3 to 0.85

For brevity, let's make a 6 X 6 X 6 steppar log (in order to do confidence limits)

XSPEC>log steppar_160059_6by6by6.log
Logging to file: steppar_160059_6by6by6.log
XSPEC>steppar 1 -0.75 1.25 6 2 0 1.0 6 3 0.3 0.85 6
 C-statistic        Delta   Mod param     Mod param     Mod param
              C-statistic         1            2            3
     221.50        47.215 -0.750000       0.00000      0.300000
...

     217.33        43.038  0.916667       1.00000      0.850000
     223.54        49.257   1.25000       1.00000      0.850000
XSPEC>log none  		!this keeps the log from getting more stuff added to it.

OK, have the .log output file, now need to modify it, and run it through Kevin's code: 

Delete the header and footer

[cohen@ruby kevin_idl]$ cp /home/kgrizza1/xspec/idl/confidence_windprof_large_show3.pro .

Load into IDLDE, compile and run; read in the log file, answer questions sent to terminal 

*the "what statistic should be plotted" question - 0,1 are slices and 2,3 are projections

Note, if you do a projection option, it'll still ask which slice of the third parameter; you have to answer, but it doesn't matter what answer you give. 

Program automatically writes out confidence_windprof_contour.ps

Making a bigger grid: 

running on amundsen - 

reran the fit, got very slightly different results: 

 ---------------------------------------------------------------------------
  Model:  windprof<1> + powerlaw<2>
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   windprof   q                 -0.297154     +/-  0.348431
    2    2    1   windprof   taustar            0.268932     +/-  0.223541
    3    3    1   windprof   umax               0.660077     +/-  0.750872E-01
    4    4    1   windprof   velocity km/s       1860.00     frozen
    5    5    1   windprof   beta                1.00000     frozen
    6    6    1   windprof   waveleng A          16.0059     frozen
    7    7    1   windprof   verbose             0.00000     frozen
    8    8    1   windprof   norm               1.546639E-04 +/-  0.286195E-04
    9    9    2   powerlaw   PhoIndex            2.00000     frozen
   10   10    2   powerlaw   norm               3.244545E-03 +/-  0.619589E-03
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 C-statistic =      174.2848     using   158 PHA bins.

XSPEC>goodness 1000
 56.10% of realizations have a fit statistic <      174.3

Doing a 30x30x30 steppar run (q on -0.8:2.0, tstar on 0.:1.0, umax on 0.25:0.9)

XSPEC>log steppar_160059_30by30by30.log
Logging to file: steppar_160059_30by30by30.log
XSPEC>steppar 1 -0.8 2.0 30 2 0. 1.0 30 3 0.25 0.9 30

running at 12:58 PM...

On Sunday, realized (after plotting above results) that we need a bigger q-grid, so, made: 

Steppar_160059_30by30by30new.log
OK…still hits upper q limit (3.5), but we’ll call that unconstrained (and umax isn’t changing as q gets higher; see contour plot)

-------------------
Plotting the spectrum and model using Kevin's programs - see dir kevin_idl/

[cohen@ruby kevin_idl]$ ls
kevin_fitplot_xspec.pro  psclose.pro  psopen.pro  windprof_160059.qdp

break the qdp file into two .dat files

[cohen@ruby kevin_idl]$ cp windprof_160059.qdp windprof_160059_p1.dat
[cohen@ruby kevin_idl]$ cp windprof_160059.qdp windprof_160059_p2.dat
[cohen@ruby kevin_idl]$ mv windprof_160059_p2.dat windprof_160059_m1.dat

edit these two files, deleting the three-line header and the unused half of the data 
NO blank lines at the end either

Running from within idlde:

open kevin_fitplot_xspec 

Some things are hardwired - xrange choice is automatic; output ps file name, etc. 
-----------------