
MLCS2k2 observer-frame fitter   

Version 0.07  Tue May  5 23:28:10 EDT 2009  sj


to run you need to set the environment variable MLCS2K2_BASEDIR to
the directory where the file you are reading is located, e.g.

% setenv MLCS2K2_BASEDIR /home/saurabh/sn/mlcs2k2.v007


all the major routines are in the pro/ subdirectory. you should either
put this near the beginning of your IDL path, or else just run IDL from
inside the pro/ subdirectory, e.g.

% cd $MLCS2K2_BASEDIR/pro
% idl
IDL>

You will also need the IDL ASTRO routines somewhere in your IDL path
(so things like mrdfits.pro, etc.). These can be downloaded from
http://idlastro.gsfc.nasa.gov/

The code has been tested on IDL 6.2, 6.3, and 7.0; let me know if you
have problems running it.

observer frame light curves go in data/
put all the useful data into one file called <SN-name>.dat
the columns should be passband, hjd, mag, err

look at the examples there to see the format

your passband names should match the available passbands (see the
.fits files located in aux/passbands). if you have a passband that's
not there you will need to create its .fits file. ask saurabh for now.


the aux/ directory also contains other details that you don't
need to worry about right now. the directories under aux/spectra
are the names of the spectral samples we have defined. 


once you have your light curve (with known passbands) in the data/
directory, you need to add its entry to the sn.info file.
this is located in the fit/ subdirectory.

add a line for your SN -- the required information is:
SN name, heliocentric redshift, *guess* at T_Bmax, 
Galactic E(B-V) from SFD maps, spectral sample to use for
K-corrections (one of the samples in aux/spectra/), prior
procedure name, and vector name. 

The guess at T_Bmax doesn't have to be great because we iterate
on it by default. Within 5 or 7 rest-frame days is good.

The chisq program calls a procedure called "prior_<name>" located 
in the fit/priors/ directory. We will add other choices for
the priors later, and we'll use this to include other priors
like rest-frame IR data, snapshot spectral information, early-time
non-detections, etc. Other possibilities (like a flat A_V, etc.)
exist in that directory now; look and see if there's something
you like.

The vecname should be one of the vectors listed in the train/
subdirectory. "final" matches the low-z rest-frame fitter
(with proper training); "early" extends "final" to -20 days.
"tweaked" adjusts the colors a bit to yield proper color
zeropoints. These vectors come with various choices of the
associated S matrix (model error matrix). "sorig" is the
raw S matrix from the training; "slowz" is the same as
"sorig" but scaled up by ~2 to better match the dispersion
of the low-z objects (and what was used in the low-z 
rest-frame fitter). Others are available (see train/README);
recommended vectors are "early-smix" or "early-slowz" or
"rv19-early-smix"

okay, that should be it. start up IDL (from the pro/ directory
if it's not in your path), and run

IDL> fit,'snname'

e.g.,  IDL> fit,'2002fw'

and the program will go to work. 

if the K-correction and Gal-ext lookup tables are not present
for that redshift/color/filter combination, they will be created.
This can be a little slow, but it's only done the first time you fit
the object.

Then we do an amoeba fit with all parameters free.

Then we do 3 iterations of a fit on a grid, with the grid points and
separations updated in each iteration to better sample the
probability. If the fit host-extinction A_V is less than 0.5 mag, we
don't fit for R_V because it is pointless (we just set R_V = 3.1).
These numbers can be changed using the "minav" and "rvminav"
parameters in defaultsrvfit.pro. You get even more flexibility if
you fix the grid, see below under the advanced usage section.

The results of the fit are in the fit/out/snname/ directory.  There
are a bunch of files. There are three log files (.fit.log, .mcs.log
and .info.log) that summarize the results and the minimum chisq.
Then there are files for each iteration; the sav files have all the
IDL information, the .fit.out files show the results (same as in the
.fit.log file). 

The .obsframe.dat files gives the model light curves in all of the
observer-frame passbands.  It also gives the model uncertainty
(S-matrix diagonal elements, plus the K-correction and Galactic
extinction uncertainties).  We also explicitly list the K-correction
used and the amount of Galactic extinction.  

The .restframe.dat files have the rest-frame magnitudes, i.e.  the
templates with that delta and corrected for host-galaxy extinction.
No K-correction or Milky Way extinction is applied.  The numbers are
given in absolute magnitudes (for H0 ~ 65; the distance modulus is not
added) sampled on the rest-frame phases.

The .datamodel.dat files show the data and the observer-frame model (now
interpolated to the times of the data), along with residuals,
diagonal elements of the C and S matrices, K-corrections and MW
extinction corrections that were applied. 

The .datamodel.dat file also calculates the SNR of each data point
(obs), as well as what the SNR of that observation *would have been*
if the noise level was the same, but the signal level was equal to
the model prediction. So for example, if a measurement is 25.0 +/-
0.1 mag, the obs SNR is about 10, but if the model predicts that the
SN was really mag 26 at that time, the model SNR is only 4, and it's
possible that the observation was a randomly high fluctuation and
using the observed signal gives a spuriously high SNR. This allows
you to set a threshold SNR (in defaultssnr.pro, currently 4) below
which to flag photometry points, using the obs SNR or the model SNR.
A "cleaned" version of the data file is output in .hisnr.dat, with
points with an obs or model SNR too low commented out. The fitter
DOES NOT use this file by default, YOU have to replace/edit the
light curve in data/

The .hisnr.dat file also excludes 5-sigma outliers (also set
in defaultssnr.pro). Again, the fitter DOES NOT LOOK AT
the .hisnr.dat file. If you want to use it, you have to
replace the light curve in data/  

Finally, the postscript files show plots of the fits (all
iterations) and the 1-d and 2-d probability distributions (grid
iterations only).

Note that the distance moduli are only *approximately* on an 
H0 ~ 65 scale. You should really treat these as relative distances only.

That's it -- basically the results are in the last line of
the .fit.log file. If things haven't converged, you can
just keep running 
IDL> fit, 'snname' 
(perhaps also with the niter keyword) and it will keep doing 
another iteration.



advanced usage


You can use the "sn.grid" file in the fit/ directory to specify
the parameters for the grid. The format of that file is 
SN name, followed by 3 numbers for each parameter. The numbers
are: # of grid points, center of the grid, distance between points

so for example, the entry:

#  SN            #mu0  mu0  mu0res     #av0  av0  av0res     #del  del  delres     #t0  t0  t0res     #rv  rv  rvres
sn99cl.flwoapo    11  30.59  0.10       11   2.5   0.10       7   0.05  0.12        5 51342.3 0.4      9  2.25  0.15

means that we do a 5-dimensional grid, with 11 mu0 points,
centered at 30.59 separated by 0.10. So that means the mu0
array goes from 30.09 to 31.09 in steps of 0.10. Similarly
for the other parameters. It's good to use an odd number
of grid points, so that the center value is actually in
the grid.

The total size of the grid is just the product of number of 
points in each direction, e.g. 11*11*7*5*9 = 38115 points.
In the current slow version of the code, you should
keep the total grid < 1 million points otherwise it will
just be too slow. < 100,000 is probably even more sane.

There are two special values of the number of grid points:

if the number of grid points is listed as 0, the program
reverts back the the default grid for that parameter (i.e.
the same as if there were no entry in the grid file for 
that SN). The grid is then automatically updated with each
iteration. Most of the time, that is what you want. 

if the number of grid points is listed as 1, the program
fixes that parameter to the center value given (even for
the amoeba). So for example,

#  SN            #mu0  mu0  mu0res     #av0  av0  av0res     #del  del  delres     #t0  t0  t0res     #rv  rv  rvres
sn98bu             0    0     0          0    0     0         0     0     0         0    0    0        1  3.10  0.00

fixes R_V for 1998bu to 3.1, while all the other paramters
are fit and their grids are automatically updated. You could
use this to fix A_V=0, or R_V=2.7, etc.

Again, you don't -need- an entry in the sn.grid file for
your SN; only if the automatic gridding isn't doing what you
want.





By default the plots just include the observer-frame passbands
which had actual data points (from the .dat file in the 
data/ directory). If you want other passbands to be plotted
(say you just want to see what the U-band model light curve
looks like, but you don't have any U-band data), you can
run:

IDL> fit, 'snname', addobspb=['U']

This will only work if the observer frame passbands you want 
correspond to rest-frame U,B,V,R or I.



Sometimes the amoeba falls into a local minimum. Usually
it's close enough to a true minimum that the subsequent grid
iterations find robust means and errors. If you are paranoid,
you can use the results of the last grid iteration to serve
as the starting point for a new amoeba. To do that:

IDL> fit, 'snname', /restart

This will then do a new amoeba (starting at the last mean
fit position), and then 3 grids. If you want just a new
amoeba and say, 2 grids, you can do:

IDL> fit, 'snname', /restart, niter=2


There is a routine called genlc to create fake light curves,
given parameters and desired observer passbands, and 
(optionally) dates or phases.

