Skip to content

Fitting pair separations across stars as a function of stellar parameters

Jump to bottom
Daniel Berke edited this page Aug 24, 2022 · 4 revisions

This page explains the process of measuring pair separations in stars, and fitting them across stars as a function of stellar effective temperature, metallicity, and surface gravity. To start with, we need a brief introduction to the star of the show, the Star class.

Overview of varconlib.Star

The Star class is used to hold all the information about a star and the fits performed on observations of it. Any information that pertains to all transitions or pairs across all observations is stored in NumPy arrays, making operations like applying CCD-position dependent wavelength corrections to measured wavelengths as simple as subtracting two arrays. Each row corresponds to data from a single observation, while each columns contains data from a single transition (or pair) across all observations. (Individual instances of pairs/transitions fitted twice in adjacent echelle orders have their own columns.) The number of columns in the array is thus the number of transitions or pairs (depending on the array) which are found in the final_transitions_selection.pkl and final_pairs_selection.pkl files, while the number of rows is equal to the number of observations of the star.

To help with finding the data for specific transitions or pairs, Star contains the following useful methods, each of which come in a "transition" and "pair" version:

  • t_label, p_label: take an integer and returns the label of the transition/pair at that row of an array. E.g. (where hd117618 is an instance of Star constructed using data from HD 117618): 
    $ hd117618.t_label(0)
>>> '4205.650Cr1_15'
$ hd117618.p_label(0)
>>> '4205.650Cr1_4210.786Fe1_15'
    Note that the number of transitions and pairs are not equal for HARPS data (and not likely to be in general). The lengths can be found as: 
    $ len(hd117618.transitionsList)
>>> 295
$ len(hd117618.pairsList)
>>> 651
  • t_index, p_index: take a transition or pair label, and return the index number of the column associated with it. E.g.: 
    $ hd117618.t_index('4937.713Cr1_37')
>>> 199
$ hd117618.p_index('4683.218Ti1_4684.870Fe1_30')
>>> 165

Initializing a Star

  1. Measure the wavelength of each transition in all observations using find_transitions.py, which uses the expected wavelengths from BERV + stellar RV to find the absorption features to measure. Fit results are saved as pickled (and compressed) GaussianFit objects.
  2. When an instance of Star is created by compiling the results for all observations from a star the following arrays are created:
    • fitMeansArray: the raw measured wavelength of each transition.
    • fitErrorsArray: the uncertainty for each fit.
    • fitOffsetsArray: the difference in velocity between the measured wavelength and the wavelength as naïvely calculated from the transition energy and the radial velocity Doppler shift.
    • ccdCorrectionsArray: the correction to be applied to each measured wavelength based on where it was measured on the CCD, from the work of Milacović et al. (2020).
    • fitOffsetsCorrectedArray: created by subtracting the Milaković corrections from the velocity offsets, fitOffsetsArray - ccdCorrectionArray.
    • fitMeansCCDCorrectedArray: created by subtracting the Milaković corrections from the raw measured wavelengths, fitMeansArray - ccdCorrectionArray.
    • obsRVOffsetsArray: a 1D array containing the median of the values of each row in fitOffsetsCorrectedArray. What this does is find a normalization factor for each observation that is sensitive to changes in radial velocity of the star between observation, such as from exoplanets.
    • fitOffsetsNormalizedArray: created by subtracted the normalization factor for each observation from all transitions offsets in that observation, e.g., fitOffsetsCorrectedArray - obsRVOffsetsArray. This removes changes in radial velocity from transition offsets. This step is not applied to the raw wavelengths, since the differential pair separations are not affected by changes in radial velocity.

Measuring pair separations

At this point it is assumed that you have run find_transitions.py on a number of stars (which may have differing numbers of observations), and have the results for each star saved as compressed (using LZMA compression) pickle file. This process will mainly use three scripts, multi_fit_stars.py (to perform fitting of transitions and pairs across stars), modify_stars.py (a convenience wrapper around methods in the Star class), and generate_stellar_database.py (which collates information across stars into one file/data structure for easy access).

  1. The first step is to run generate_stellar_database.py to collect data from stars together in once place for multi_fit_stars.py to read. The basic call will look like this: 
    $ generate_stellar_database.py <main_directory> <names_of_star_directories_to_use> --transitions-only --casagrande2011
    Here, <main_directory> should be the directory containing the directories for individual stars, likely the same as that defined by data_output_dir in the VarConLib config file. --transitions-only is because data on pair separations doesn't exist yet, so trying to read it from Star will cause a crash. --casagrande2011 refers to the paper to use for stellar parameter values (the other option is --nordstrom2004, but it should be avoided). An example might look like: 
    $ generate_stellar_database.py /Users/dberke/data_output Vesta HD45184 HD146233 HD20782 --transitions-only --casagrande2011
    If a Star instance has not been created for each star yet (and its data saved to HDF5) they will be created now. There is also a --recreate-stars flag which may be used to force the rebuilding of all named stars. The script will save its output to the directory defined as databases_dir in the VarConLib config as "stellar_db_uncorrected.hdf5". (Despite the usage of the term 'database', there is no SQL involved.)
  2. Next, run multi_fit_stars.py, which will use the data collected in the file from generate_stellar_database.py. multi_fit_stars.py takes several important flags:
    • (Required) A flag denoting whether it is to fit transitions or pairs. This can be either -T, --transitions, or -P, --pairs. It is necessary to have one of these flags for it to run, but not both.
    • (Required) A flag denoting the form of the function used to fit the transitions/pairs. All functions are multi-variate in stellar effective temperature, metallicity, and surface gravity. For our purposes, we found a function quadratic in all terms (--quadratic) to work best.
    • (Optional) A flag determining the sigma-clipping cutoff limit. The default is 2.5. For this step, use the following call:
    $ multi_fit_stars.py -T --quadratic --sigma 3.0
    This script will, for each transition in final_transitions_selection.pkl, take the weighted mean value for each star (created from all observations of that star) and then fit these weighted mean values as a function of each star's three stellar parameters of interest. The fitting results, with the function used and the best-fit coefficients for each transition, will be saved to an HDF5 file in a directory called fit_params within the data_output_dir defined in the VarConLib config. The name of the file will be of the form "<model_name><fit_target><#sigma>sigma_hdf5", e.g., "quadratic_transitions_3.0sigma.hdf5". There will also be a directory created within the data_output_dir with a name of the form stellar_parameter_fits_<fit_target>_<#sigma>sigma/<model_name>, e.g. stellar_parameter_fits_transitions_3.0sigma/quadratic. This folder will contain plots of the fit for each transition, and a CSV file containing some information about the fits for each transition, namely the chi-squared value, sigma, and sigma_** for both pre- and post- observations.
  3. After this, there are three methods of Star which will need to be called (though not all at once): createTransitionModelCorrectedArrays, createPairSeparationArrays, and createPairModelCorrectedArrays. For now, the first two methods need to be called in each star, which can be accomplished with the modify_stars.py script. This script takes a list of star names (assumed to correspond to directories in the data_output_dir and a few flags. For this step, use the following parameters: 
    $ modify_stars.py <list of star names> --transitions --sigma 3.0
    This will go through each star, create a Star object, and run createTransitionModelCorrectedArrays and createPairSeparationArrays. By default, this will use a function quadratic in effective temperature, metallicity, and surface gravity. (An additional flag, --linear, can be provided to use a purely linear function.) The first method will create the following arrays and add them to the saved information for that Star:
    • transitionModelOffsetsArray: contains the offset from the model for each transition measurement for the star, for all observations. Any values which are found to be outliers beyond the 'sigma' value given will be listed as 'NaN', so that they cannot be accidentally used in future. Generally, unless stated otherwise, values are in m/s; all arrays should also have units attached using the "Unyt" Python package.
    • transitionModelErrorsArray: a copy of fitErrorsArray, masked so that the same values as transitionModelOffsetsArray are NaN.
    • transitionModelArray: this array has only two rows, for pre- (0) and post (1) fiber-change observations, with as many columns as transitions. Each column contains the model value(s) for the associated transition.
    • transitionSysErrorsArray: similarly, this array has 2 rows and as many columns as transitions, and contains the sigma_** value(s) for each transition for the pre- and post- eras. After this, createPairSeparationArrays will create the following arrays (this is where the actual pair separations are measured):
    • pairSeparationsArray: contains the measured pair separation for each pair in final_pair_selection.pkl, for each observation of the star. These are measured from the values in firMeansCCDCorrectedArray (so the normalization step to remove changes in stellar radial velocity from e.g. exoplanets is not applied to these values). During this step fitMeansCCDCorrectedArray has the mask from createTransitionModelCorrectedArrays applied, so any outliers are marked as NaN.
    • pairSepErrorsArray: the quadrature sum of the error for both component transitions in the pair, from fitErrorsArray. At the end of this step the pair separation measurements will have been taken, and are ready for fitting similar to the transitions.
  4. At this stage, run generate_stellar_database.py again without the --transitions-only flag, e.g.: 
    $ generate_stellar_database.py /Users/dberke/data_output Vesta HD45184 HD146233 HD20782 --casagrande2011
    This will collect the pair separation measurements into the database HDF5 file (overwriting the previous one).
  5. Now, multi_fit_stars.py can be run again to fit pair separations across all stars. For pairs, we decided on a slightly more conservative 4-sigma cut: 
    $ multi_fit_stars.py -P --quadratic --sigma 4.0
    The script will produce output similar to that in step 2 from fitting transitions, just with 'pairs' in the name rather than 'transitions'.
  6. Finally, run modify_stars.py again with the following parameters: 
    $ modify_stars.py <list of star names> --pairs --sigma 4.0
    This will call the createPairModelCorrectedArrays method of each Star, creating some arrays analogous to those produced for createTransitionModelCorrectedArrays:
    • pairModelOffsetsArray: an array containing the measured pair separation minus the model value for each pair in the star. Any pair separations which were not created due to one or both transitions being outliers, or which are themselves outliers, are denoted by NaN.
    • pairModelErrorsArray: a copy of pairSepErrorsArray, masked the same as pairModelOffsetsArray.
    • pairModelArray: similar to transitionModelArray, this is a 2-row table containing the model-expected value for each pair, for the pre- and post- fiber-change eras. (Pre = 0, post = 1.)
    • pairSysErrorsArray: similar to transitionSysErrorsArray, this is a 2-row table containing the sigma_** (excess scatter) for each pair.
  7. At this point, all the information is stored in Star objects (if I'd thought of it and had the time, I'd probably have made it so that running generate_stellar_database.py a third time would collect all the pair model comparison information). To write out results for stars, use the save_pair_separations.py script like so: 
    $ save_pair_separations.py <list of star names in 'data_output_dir'> --pairs

This will create a directory within data_output_dir named pair_separation_files, within which will be two directories named "pre" and "post". Inside these directories the script will create a .csv file for each pair, with one row for each star given on the command line. These rows will contain the star name, number of observations in the relevant era (pre or post), the offset from the model for weighted mean offset of that pair, the statistical and systematic error, the chi-squared per degree of freedom for that pair over all observations of the star, and then the offset, statistical error, systematic error, and chi-squared per DOF for the two transitions in the pair (blue, then red). (This file is generated by calling the formatPairData method of Star.)

And that should be it. The data is of course also all still available for individual stars down to individual measurements, allowing for very fine-grained analysis if desired. I'll add more to this if I've forgotten something.

Clone this wiki locally