Basic COMPAS data analysis¶
Most of the post-processing material presented here is written in python, and makes extensive use of the numpy package for rapid computation on large data arrays. Here we show two important basics of python/numpy in the context of investigating your COMPAS simulation.
Material¶
1) Inspecting the data¶
Look at the data to see which parameters are available and check that it matches expectations.
2) Slicing the data¶
Select specific systems and their parameters using seeds.
3) Visualizing the data¶
Binning and visualising your data.
[ ]:
For the following sections, you will need to have the following packages installed.¶
numpy, h5py, time, matplotlib¶
[3]:
#python libraries
import os, sys
import numpy as np # for handling arrays
import h5py as h5 # for reading the COMPAS data
import time # for finding computation time
import matplotlib.pyplot as plt #for plotting
# Import COMPAS specific scripts
compasRootDir = os.environ['COMPAS_ROOT_DIR']
sys.path.append(compasRootDir + 'postProcessing/PythonScripts')
from compasUtils import printCompasDetails, getEventHistory, getEventStrings
# Choose an output hdf5 file to work with
pathToData = 'COMPAS_Tutorial_Output.h5'
# This is known as an ipython magic command, and allows plots to be produced within the notebook
%matplotlib inline
[ ]:
1. Inspecting the data¶
Often the first thing you want to do with new data is simply to look at it! Getting familiar with the data, including available parameters, size of the data file, etc. will help to inform how best to proceed with the analysis. We provide several useful functions for inspecting the data, printCompasDetails, getEventHistory, and getEventStrings.
–
If you want to create an alternative COMPAS_Output.h5, see Section 1 Working With HDF5, or download some data from our Zenodo database.
Note: These cells may take a long time if you test them on large datasets.
[4]:
Data = h5.File(pathToData)
print(list(Data.keys()))
['BSE_Common_Envelopes', 'BSE_Double_Compact_Objects', 'BSE_RLOF', 'BSE_Supernovae', 'BSE_System_Parameters', 'Run_Details']
The output above represents the event categories available from the particular run. If you used the output produced in the previous tutorial, you should see ['BSE_Common_Envelopes', 'BSE_Double_Compact_Objects', 'BSE_RLOF', 'BSE_Supernovae', 'BSE_System_Parameters', 'Run_Details']. Note that for smaller runs which do not produce any of a particular type of output, the output category will not be created.
Brief description of the categories: - ‘BSE_System_Parameters’: Initial state of the binary - ‘BSE_RLOF’: Any mass transfer events that occured within the binary - ‘BSE_Common_Envelopes’: If any of the mass transfer events were unstable, details will be included here. - ‘BSE_Supernovae’: Parameters and outcome of any supernovae that occured in the binary - ‘BSE_Double_Compact_Objects’: Includes key information of all binaries which end their lives as an intact pair of compact obects (either neutron stars or black holes) - ‘Run_Details’: Information on the input settings supplied to the Compas run
To extract the data from these categories, we use the following syntax
[5]:
SPs = Data['BSE_System_Parameters']
MTs = Data['BSE_RLOF']
CEs = Data['BSE_Common_Envelopes']
SNe = Data['BSE_Supernovae']
DCs = Data['BSE_Double_Compact_Objects']
Each of these is a dictionary mapping parameter names (keys) to an array of values
[6]:
print(SPs.keys())
<KeysViewHDF5 ['Add_Options_To_SysParms', 'Allow_MS_To_Survive_CE', 'Allow_RLOF@Birth', 'Allow_Touching@Birth', 'BB_Mass_xFer_Stblty_Prscrptn', 'BH_Kicks', 'CE_Alpha', 'CE_Alpha_Thermal', 'CE_Lambda', 'CE_Lambda_Multiplier', 'CE_Lambda_Prscrptn', 'CE_Mass_Accr_Constant', 'CE_Mass_Accr_Max', 'CE_Mass_Accr_Min', 'CE_Mass_Accr_Prscrptn', 'CE_Recomb_Enrgy_Dnsty', 'CE_Slope_Kruckow', 'CHE_Mode', 'CH_on_MS(1)', 'CH_on_MS(2)', 'Check_Photon_Tiring_Limit', 'Circularise@MT', 'Conserve_AngMom@Circ', 'Cool_WindMassLoss_Multipl', 'Eccentricity', 'Eccentricity@ZAMS', 'Eccentricity_Dstrbtn', 'Eccentricity_Dstrbtn_Max', 'Eccentricity_Dstrbtn_Min', 'Eddington_Accr_Factor', 'Envelope_State_Prscrptn', 'Equilibrated_At_Birth', 'Error', 'Evolution_Mode', 'Fryer_SN_Engine', 'Initial_Mass', 'Initial_Mass(1)', 'Initial_Mass(2)', 'Initial_Mass_Func_Max', 'Initial_Mass_Func_Min', 'Initial_Mass_Func_Power', 'Initial_Mass_Function', 'Kick_Direction_Dstrbtn', 'Kick_Direction_Power', 'Kick_Magnitude', 'Kick_Magnitude(1)', 'Kick_Magnitude(2)', 'Kick_Magnitude_Dstrbtn', 'Kick_Magnitude_Dstrbtn_Max', 'Kick_Magnitude_Random', 'Kick_Magnitude_Random(1)', 'Kick_Magnitude_Random(2)', 'Kick_Mean_Anomaly(1)', 'Kick_Mean_Anomaly(2)', 'Kick_Phi(1)', 'Kick_Phi(2)', 'Kick_Scaling_Factor', 'Kick_Theta(1)', 'Kick_Theta(2)', 'LBV_Factor', 'LBV_Mass_Loss_Prscrptn', 'MCBUR1', 'MM_Kick_Multiplier_BH', 'MM_Kick_Multiplier_NS', 'MT_Acc_Efficiency_Prscrptn', 'MT_AngMom_Loss_Prscrptn', 'MT_Fraction_Accreted', 'MT_JLoss', 'MT_Rejuvenation_Prscrptn', 'MT_Thermal_Limit_C', 'MT_Thermally_Lmtd_Variation', 'Mass@ZAMS(1)', 'Mass@ZAMS(2)', 'Mass_Loss_Prscrptn', 'Mass_Ratio', 'Mass_Ratio_Dstrbtn', 'Mass_Ratio_Dstrbtn_Max', 'Mass_Ratio_Dstrbtn_Min', 'Max_Donor_Mass', 'Max_Evolution_Time', 'Max_NS_Mass', 'Max_Timesteps', 'Merger', 'Merger_At_Birth', 'Metallicity', 'Metallicity@ZAMS(1)', 'Metallicity@ZAMS(2)', 'Metallicity_Dstrbtn', 'Metallicity_Dstrbtn_Max', 'Metallicity_Dstrbtn_Min', 'Min_Secondary_Mass', 'NS_EOS', 'Neutrino_Mass_Loss_Assmptn', 'Neutrino_Mass_Loss_Value', 'Omega@ZAMS(1)', 'Omega@ZAMS(2)', 'Orbital_Period', 'Orbital_Period_Dstrbtn', 'Orbital_Period_Max', 'Orbital_Period_Min', 'Overall_WindMassLoss_Multipl', 'PISN_Lower_Limit', 'PISN_Upper_Limit', 'PPI_Lower_Limit', 'PPI_Prscrptn', 'PPI_Upper_Limit', 'Pulsar_Mag_Field_Decay_mScale', 'Pulsar_Mag_Field_Decay_tScale', 'Pulsar_Mag_Field_Dstrbtn', 'Pulsar_Mag_Field_Dstrbtn_Max', 'Pulsar_Mag_Field_Dstrbtn_Min', 'Pulsar_Minimum_Mag_Field', 'Pulsar_Spin_Period_Dstrbtn', 'Pulsar_Spin_Period_Dstrbtn_Max', 'Pulsar_Spin_Period_Dstrbtn_Min', 'Remnant_Mass_Prscrptn', 'Rotational_Frequency', 'Rotational_Frequency(1)', 'Rotational_Frequency(2)', 'Rotational_Velocity_Dstrbtn', 'SEED', 'SEED(CMDLINE)', 'SEED(OPTION)', 'SN_Kick_Magnitude_Random_Number(1)', 'SN_Kick_Magnitude_Random_Number(2)', 'Semi-Major_Axis', 'Semi-Major_Axis_Dstrbtn', 'Semi-Major_Axis_Dstrbtn_Max', 'Semi-Major_Axis_Dstrbtn_Min', 'Semi-Major_Axis_Dstrbtn_Power', 'SemiMajorAxis@ZAMS', 'Sigma_Kick_CCSN_BH', 'Sigma_Kick_CCSN_NS', 'Sigma_Kick_ECSN', 'Sigma_Kick_USSN', 'Stellar_Type(1)', 'Stellar_Type(2)', 'Stellar_Type@ZAMS(1)', 'Stellar_Type@ZAMS(2)', 'Stellar_Zeta_Prscrptn', 'Unbound', 'WR_Factor', 'Zeta_Adiabatic_Arbitrary', 'Zeta_Main_Sequence', 'Zeta_Radiative_Envelope_Giant']>
One of the most important parameters in the COMPAS output is the system seed. The seed represents the unique identifier to a specific system in a simulation. It is also used as the seed value in random number generation, which is useful when trying to reproduce a given system identically.
If we want to view, say, the random seeds in the system parameters file, we run
[7]:
seedsSP = SPs['SEED'][()]
print(seedsSP)
[1636090318 1636090389 1636091116 101048 102660]
printCompasDetails¶
This is useful for extracting the arrays of single parameters, but for a more convenient view of the whole system parameters file, we can use the printCompasDetails function.
[8]:
printCompasDetails(SPs) # Note - the output of this is a pandas dataframe
[8]:
| SEED | (units) | 1636090318 | 1636090389 | 1636091116 | 101048 | 102660 |
|---|---|---|---|---|---|---|
| Add_Options_To_SysParms | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Allow_MS_To_Survive_CE | Flag | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Allow_RLOF@Birth | Flag | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Allow_Touching@Birth | Flag | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| BB_Mass_xFer_Stblty_Prscrptn | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| BH_Kicks | - | 3.000000e+00 | 3.000000e+00 | 3.000000e+00 | 3.000000e+00 | 3.000000e+00 |
| CE_Alpha | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| CE_Alpha_Thermal | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| CE_Lambda | - | 1.000000e-01 | 1.000000e-01 | 1.000000e-01 | 1.000000e-01 | 1.000000e-01 |
| CE_Lambda_Multiplier | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| CE_Lambda_Prscrptn | - | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 |
| CE_Mass_Accr_Constant | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| CE_Mass_Accr_Max | Msol | 1.000000e-01 | 1.000000e-01 | 1.000000e-01 | 1.000000e-01 | 1.000000e-01 |
| CE_Mass_Accr_Min | Msol | 4.000000e-02 | 4.000000e-02 | 4.000000e-02 | 4.000000e-02 | 4.000000e-02 |
| CE_Mass_Accr_Prscrptn | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| CE_Recomb_Enrgy_Dnsty | erg g^-1 | 1.500000e+13 | 1.500000e+13 | 1.500000e+13 | 1.500000e+13 | 1.500000e+13 |
| CE_Slope_Kruckow | - | -8.333333e-01 | -8.333333e-01 | -8.333333e-01 | -8.333333e-01 | -8.333333e-01 |
| CHE_Mode | - | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 |
| CH_on_MS(1) | State | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| CH_on_MS(2) | State | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Check_Photon_Tiring_Limit | Flag | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Circularise@MT | Flag | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Conserve_AngMom@Circ | Flag | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Cool_WindMassLoss_Multipl | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Eccentricity | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Eccentricity@ZAMS | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Eccentricity_Dstrbtn | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Eccentricity_Dstrbtn_Max | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Eccentricity_Dstrbtn_Min | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Eddington_Accr_Factor | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Envelope_State_Prscrptn | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Equilibrated_At_Birth | Event | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Error | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Evolution_Mode | Mode | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Fryer_SN_Engine | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Initial_Mass | Msol | 5.000000e+00 | 5.000000e+00 | 5.000000e+00 | 5.000000e+00 | 5.000000e+00 |
| Initial_Mass(1) | Msol | 1.270356e+01 | 7.655521e+01 | 1.134459e+02 | 2.494177e+01 | 2.674675e+01 |
| Initial_Mass(2) | Msol | 2.403020e+00 | 3.954175e+01 | 8.169763e+01 | 1.542647e+01 | 1.657629e+01 |
| Initial_Mass_Func_Max | Msol | 1.500000e+02 | 1.500000e+02 | 1.500000e+02 | 1.500000e+02 | 1.500000e+02 |
| Initial_Mass_Func_Min | Msol | 5.000000e+00 | 5.000000e+00 | 5.000000e+00 | 5.000000e+00 | 5.000000e+00 |
| Initial_Mass_Func_Power | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Initial_Mass_Function | - | 3.000000e+00 | 3.000000e+00 | 3.000000e+00 | 3.000000e+00 | 3.000000e+00 |
| Kick_Direction_Dstrbtn | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Kick_Direction_Power | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Kick_Magnitude | kms^-1 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Kick_Magnitude(1) | kms^-1 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Kick_Magnitude(2) | kms^-1 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Kick_Magnitude_Dstrbtn | - | 3.000000e+00 | 3.000000e+00 | 3.000000e+00 | 3.000000e+00 | 3.000000e+00 |
| Kick_Magnitude_Dstrbtn_Max | - | -1.000000e+00 | -1.000000e+00 | -1.000000e+00 | -1.000000e+00 | -1.000000e+00 |
| Kick_Magnitude_Random | - | 7.463209e-01 | 2.867494e-01 | 5.613971e-01 | 6.265861e-01 | 9.722633e-01 |
| Kick_Magnitude_Random(1) | - | 7.690129e-01 | 1.203533e-01 | 5.599079e-01 | 1.391730e-02 | 1.704616e-01 |
| Kick_Magnitude_Random(2) | - | 3.773938e-01 | 5.905514e-01 | 8.641669e-01 | 8.270186e-01 | 1.328111e-01 |
| Kick_Mean_Anomaly(1) | - | 5.189746e+00 | 1.824014e+00 | 9.406175e-01 | 5.294177e+00 | 3.100337e+00 |
| Kick_Mean_Anomaly(2) | - | 2.995799e+00 | 5.662449e+00 | 1.977927e+00 | 2.189489e+00 | 5.305596e+00 |
| Kick_Phi(1) | - | 6.485718e-01 | -6.633342e-01 | -1.781796e-01 | -7.482819e-01 | 1.205367e+00 |
| Kick_Phi(2) | - | 2.182405e-02 | 3.901919e-01 | 5.983047e-01 | 9.572699e-02 | -3.405070e-01 |
| Kick_Scaling_Factor | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Kick_Theta(1) | - | 3.478068e+00 | 3.615591e+00 | 1.748962e+00 | 5.110205e+00 | 5.544482e+00 |
| Kick_Theta(2) | - | 2.636540e-01 | 5.426292e-01 | 9.971980e-01 | 5.717921e+00 | 3.657304e+00 |
| LBV_Factor | - | 1.500000e+00 | 1.500000e+00 | 1.500000e+00 | 1.500000e+00 | 1.500000e+00 |
| LBV_Mass_Loss_Prscrptn | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| MCBUR1 | Msol | 1.600000e+00 | 1.600000e+00 | 1.600000e+00 | 1.600000e+00 | 1.600000e+00 |
| MM_Kick_Multiplier_BH | - | 2.000000e+02 | 2.000000e+02 | 2.000000e+02 | 2.000000e+02 | 2.000000e+02 |
| MM_Kick_Multiplier_NS | - | 4.000000e+02 | 4.000000e+02 | 4.000000e+02 | 4.000000e+02 | 4.000000e+02 |
| MT_Acc_Efficiency_Prscrptn | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| MT_AngMom_Loss_Prscrptn | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| MT_Fraction_Accreted | - | 5.000000e-01 | 5.000000e-01 | 5.000000e-01 | 5.000000e-01 | 5.000000e-01 |
| MT_JLoss | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| MT_Rejuvenation_Prscrptn | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| MT_Thermal_Limit_C | - | 1.000000e+01 | 1.000000e+01 | 1.000000e+01 | 1.000000e+01 | 1.000000e+01 |
| MT_Thermally_Lmtd_Variation | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Mass@ZAMS(1) | Msol | 1.270356e+01 | 7.655521e+01 | 1.134459e+02 | 2.494177e+01 | 2.674675e+01 |
| Mass@ZAMS(2) | Msol | 2.403020e+00 | 3.954175e+01 | 8.169763e+01 | 1.542647e+01 | 1.657629e+01 |
| Mass_Loss_Prscrptn | - | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 |
| Mass_Ratio | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Mass_Ratio_Dstrbtn | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Mass_Ratio_Dstrbtn_Max | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Mass_Ratio_Dstrbtn_Min | - | 1.000000e-02 | 1.000000e-02 | 1.000000e-02 | 1.000000e-02 | 1.000000e-02 |
| Max_Donor_Mass | Msol | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 |
| Max_Evolution_Time | Myr | 1.370000e+04 | 1.370000e+04 | 1.370000e+04 | 1.370000e+04 | 1.370000e+04 |
| Max_NS_Mass | Msol | 2.500000e+00 | 2.500000e+00 | 2.500000e+00 | 2.500000e+00 | 2.500000e+00 |
| Max_Timesteps | Count | 9.999900e+04 | 9.999900e+04 | 9.999900e+04 | 9.999900e+04 | 9.999900e+04 |
| Merger | Event | 1.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Merger_At_Birth | Event | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Metallicity | - | 1.420000e-02 | 1.420000e-02 | 1.420000e-02 | 1.420000e-02 | 1.420000e-02 |
| Metallicity@ZAMS(1) | - | 1.420000e-02 | 1.420000e-02 | 1.420000e-02 | 1.420000e-02 | 1.420000e-02 |
| Metallicity@ZAMS(2) | - | 1.420000e-02 | 1.420000e-02 | 1.420000e-02 | 1.420000e-02 | 1.420000e-02 |
| Metallicity_Dstrbtn | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Metallicity_Dstrbtn_Max | - | 3.000000e-02 | 3.000000e-02 | 3.000000e-02 | 3.000000e-02 | 3.000000e-02 |
| Metallicity_Dstrbtn_Min | - | 1.000000e-04 | 1.000000e-04 | 1.000000e-04 | 1.000000e-04 | 1.000000e-04 |
| Min_Secondary_Mass | Msol | 1.000000e-01 | 1.000000e-01 | 1.000000e-01 | 1.000000e-01 | 1.000000e-01 |
| NS_EOS | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Neutrino_Mass_Loss_Assmptn | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Neutrino_Mass_Loss_Value | - | 1.000000e-01 | 1.000000e-01 | 1.000000e-01 | 1.000000e-01 | 1.000000e-01 |
| Omega@ZAMS(1) | Hz | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Omega@ZAMS(2) | Hz | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Orbital_Period | days | 1.000000e-01 | 1.000000e-01 | 1.000000e-01 | 1.000000e-01 | 1.000000e-01 |
| Orbital_Period_Dstrbtn | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Orbital_Period_Max | days | 1.000000e+03 | 1.000000e+03 | 1.000000e+03 | 1.000000e+03 | 1.000000e+03 |
| Orbital_Period_Min | days | 1.100000e+00 | 1.100000e+00 | 1.100000e+00 | 1.100000e+00 | 1.100000e+00 |
| Overall_WindMassLoss_Multipl | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| PISN_Lower_Limit | Msol | 6.000000e+01 | 6.000000e+01 | 6.000000e+01 | 6.000000e+01 | 6.000000e+01 |
| PISN_Upper_Limit | Msol | 1.350000e+02 | 1.350000e+02 | 1.350000e+02 | 1.350000e+02 | 1.350000e+02 |
| PPI_Lower_Limit | Msol | 3.500000e+01 | 3.500000e+01 | 3.500000e+01 | 3.500000e+01 | 3.500000e+01 |
| PPI_Prscrptn | - | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 |
| PPI_Upper_Limit | Msol | 6.000000e+01 | 6.000000e+01 | 6.000000e+01 | 6.000000e+01 | 6.000000e+01 |
| Pulsar_Mag_Field_Decay_mScale | Msol | 2.500000e-02 | 2.500000e-02 | 2.500000e-02 | 2.500000e-02 | 2.500000e-02 |
| Pulsar_Mag_Field_Decay_tScale | Myr | 1.000000e+03 | 1.000000e+03 | 1.000000e+03 | 1.000000e+03 | 1.000000e+03 |
| Pulsar_Mag_Field_Dstrbtn | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Pulsar_Mag_Field_Dstrbtn_Max | AU | 1.300000e+01 | 1.300000e+01 | 1.300000e+01 | 1.300000e+01 | 1.300000e+01 |
| Pulsar_Mag_Field_Dstrbtn_Min | AU | 1.100000e+01 | 1.100000e+01 | 1.100000e+01 | 1.100000e+01 | 1.100000e+01 |
| Pulsar_Minimum_Mag_Field | Gauss | 8.000000e+00 | 8.000000e+00 | 8.000000e+00 | 8.000000e+00 | 8.000000e+00 |
| Pulsar_Spin_Period_Dstrbtn | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Pulsar_Spin_Period_Dstrbtn_Max | AU | 1.000000e+02 | 1.000000e+02 | 1.000000e+02 | 1.000000e+02 | 1.000000e+02 |
| Pulsar_Spin_Period_Dstrbtn_Min | AU | 1.000000e+01 | 1.000000e+01 | 1.000000e+01 | 1.000000e+01 | 1.000000e+01 |
| Remnant_Mass_Prscrptn | - | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 |
| Rotational_Frequency | Hz | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Rotational_Frequency(1) | Hz | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Rotational_Frequency(2) | Hz | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Rotational_Velocity_Dstrbtn | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| SEED(CMDLINE) | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| SEED(OPTION) | - | 1.636090e+09 | 1.636090e+09 | 1.636091e+09 | 1.010480e+05 | 1.026600e+05 |
| SN_Kick_Magnitude_Random_Number(1) | - | 7.690129e-01 | 1.203533e-01 | 5.599079e-01 | 1.391730e-02 | 1.704616e-01 |
| SN_Kick_Magnitude_Random_Number(2) | - | 3.773938e-01 | 5.905514e-01 | 8.641669e-01 | 8.270186e-01 | 1.328111e-01 |
| Semi-Major_Axis | AU | 4.709733e+00 | 8.551670e-01 | 1.013902e+00 | 5.196680e-01 | 8.079884e-01 |
| Semi-Major_Axis_Dstrbtn | - | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| Semi-Major_Axis_Dstrbtn_Max | AU | 1.000000e+03 | 1.000000e+03 | 1.000000e+03 | 1.000000e+03 | 1.000000e+03 |
| Semi-Major_Axis_Dstrbtn_Min | AU | 1.000000e-02 | 1.000000e-02 | 1.000000e-02 | 1.000000e-02 | 1.000000e-02 |
| Semi-Major_Axis_Dstrbtn_Power | - | -1.000000e+00 | -1.000000e+00 | -1.000000e+00 | -1.000000e+00 | -1.000000e+00 |
| SemiMajorAxis@ZAMS | AU | 4.709733e+00 | 8.551670e-01 | 1.013902e+00 | 5.196680e-01 | 8.079884e-01 |
| Sigma_Kick_CCSN_BH | kms^-1 | 2.650000e+02 | 2.650000e+02 | 2.650000e+02 | 2.650000e+02 | 2.650000e+02 |
| Sigma_Kick_CCSN_NS | kms^-1 | 2.650000e+02 | 2.650000e+02 | 2.650000e+02 | 2.650000e+02 | 2.650000e+02 |
| Sigma_Kick_ECSN | kms^-1 | 3.000000e+01 | 3.000000e+01 | 3.000000e+01 | 3.000000e+01 | 3.000000e+01 |
| Sigma_Kick_USSN | kms^-1 | 3.000000e+01 | 3.000000e+01 | 3.000000e+01 | 3.000000e+01 | 3.000000e+01 |
| Stellar_Type(1) | - | 1.300000e+01 | 1.400000e+01 | 1.400000e+01 | 1.300000e+01 | 1.300000e+01 |
| Stellar_Type(2) | - | 1.000000e+00 | 1.300000e+01 | 1.400000e+01 | 1.300000e+01 | 1.300000e+01 |
| Stellar_Type@ZAMS(1) | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Stellar_Type@ZAMS(2) | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Stellar_Zeta_Prscrptn | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Unbound | State | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| WR_Factor | - | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 |
| Zeta_Adiabatic_Arbitrary | - | 1.000000e+04 | 1.000000e+04 | 1.000000e+04 | 1.000000e+04 | 1.000000e+04 |
| Zeta_Main_Sequence | - | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 | 2.000000e+00 |
| Zeta_Radiative_Envelope_Giant | - | 6.500000e+00 | 6.500000e+00 | 6.500000e+00 | 6.500000e+00 | 6.500000e+00 |
PrintCompasDetails optionally also takes seeds as arguments, to focus on specific systems.
[9]:
seedsMT = MTs['SEED'][()]
firstThreeUniqueSeeds = np.unique(seedsMT)[0:3]
print("Look at these seeds: ", firstThreeUniqueSeeds)
printCompasDetails(MTs, firstThreeUniqueSeeds)
Look at these seeds: [ 101048 102660 1636090318]
[9]:
| SEED | (units) | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | ... | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 101048 | 101048 | 101048 | 102660 | 102660 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| CEE>MT | State | 1.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 1.000000 | 0.000000 | 1.000000 | 0.000000 | 0.000000 | 1.000000 |
| Eccentricity<MT | - | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.743907 | 0.000000 | 0.707907 | 0.000000 | 0.000000 | 0.859918 |
| Eccentricity>MT | - | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.743907 | 0.000000 | 0.000000 | 0.000000 | 0.295882 | 0.000000 | 0.000000 |
| MT_Event_Counter | Count | 1.000000 | 2.000000 | 3.000000 | 4.000000 | 5.000000 | 6.000000 | 7.000000 | 8.000000 | 9.000000 | ... | 29.000000 | 30.000000 | 31.000000 | 32.000000 | 33.000000 | 1.000000 | 2.000000 | 3.000000 | 1.000000 | 2.000000 |
| Mass(1)<MT | Msol | 11.773666 | 3.943185 | 4.952215 | 4.957205 | 4.957273 | 4.957274 | 4.957274 | 4.957275 | 4.957275 | ... | 4.957279 | 4.957279 | 4.957279 | 4.957279 | 1.360533 | 23.939714 | 2.008319 | 2.008319 | 25.483363 | 2.155738 |
| Mass(1)>MT | Msol | 3.943185 | 4.952215 | 4.957205 | 4.957273 | 4.957274 | 4.957274 | 4.957275 | 4.957275 | 4.957275 | ... | 4.957279 | 4.957279 | 4.957279 | 1.360533 | 1.360533 | 7.053737 | 2.008319 | 2.008356 | 7.681503 | 2.155738 |
| Mass(2)<MT | Msol | 2.403020 | 2.403020 | 1.393990 | 1.389000 | 1.388932 | 1.388931 | 1.388931 | 1.388930 | 1.388930 | ... | 1.388926 | 1.388926 | 1.388926 | 1.388925 | 1.388925 | 15.378674 | 20.247173 | 4.779758 | 16.504957 | 20.690239 |
| Mass(2)>MT | Msol | 2.403020 | 1.393990 | 1.389000 | 1.388932 | 1.388931 | 1.388931 | 1.388930 | 1.388930 | 1.388930 | ... | 1.388926 | 1.388926 | 1.388925 | 1.388925 | 1.388925 | 20.692673 | 5.688526 | 1.587026 | 21.164807 | 5.868991 |
| RLOF(1)<MT | Event | 0.000000 | 1.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| RLOF(1)>MT | State | 1.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 1.000000 | 0.000000 | 0.000000 | 1.000000 | 0.000000 |
| RLOF(2)<MT | Event | 0.000000 | 0.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | ... | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| RLOF(2)>MT | State | 0.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | ... | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 0.000000 | 1.000000 | 1.000000 | 0.000000 | 1.000000 |
| Radius(1)<MT | Rsol | 569.606809 | 0.614182 | 0.500045 | 0.505819 | 0.511921 | 0.518248 | 0.524816 | 0.531644 | 0.538753 | ... | 0.803117 | 0.828516 | 0.856161 | 0.886326 | 0.000014 | 56.758287 | 0.000014 | 0.000014 | 78.043651 | 0.000014 |
| Radius(1)>MT | Rsol | 0.597875 | 0.494328 | 0.499924 | 0.505817 | 0.511921 | 0.518248 | 0.524816 | 0.531644 | 0.538753 | ... | 0.803116 | 0.828516 | 0.856161 | 0.000014 | 0.000014 | 0.754690 | 0.000014 | 0.000014 | 0.795761 | 0.000014 |
| Radius(2)<MT | Rsol | 1.698892 | 1.698897 | 1.357432 | 1.353888 | 1.353840 | 1.353839 | 1.353840 | 1.353840 | 1.353840 | ... | 1.353844 | 1.353845 | 1.353845 | 1.353845 | 1.353968 | 6.931882 | 182.606031 | 1.215352 | 7.379090 | 265.605381 |
| Radius(2)>MT | Rsol | 1.698892 | 1.357432 | 1.353887 | 1.353839 | 1.353839 | 1.353839 | 1.353839 | 1.353840 | 1.353840 | ... | 1.353844 | 1.353844 | 1.353845 | 1.353845 | 1.353968 | 7.487812 | 0.659753 | 0.000014 | 7.801333 | 0.672806 |
| SemiMajorAxis<MT | Rsol | 1079.189435 | 2.594059 | 4.880604 | 4.905801 | 4.906144 | 4.906150 | 4.906152 | 4.906153 | 4.906154 | ... | 4.906175 | 4.906176 | 4.906177 | 4.906178 | 2.936388 | 114.730674 | 999.159759 | 2.589867 | 179.269780 | 2934.242255 |
| SemiMajorAxis>MT | Rsol | 2.594059 | 4.880604 | 4.905801 | 4.906144 | 4.906150 | 4.906152 | 4.906153 | 4.906154 | 4.906156 | ... | 4.906176 | 4.906177 | 4.906178 | 2.936388 | 0.751990 | 282.770436 | 2.282538 | 2.230194 | 427.793422 | 7.677582 |
| Stellar_Type(1)<MT | 5.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | ... | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 13.000000 | 2.000000 | 13.000000 | 13.000000 | 2.000000 | 13.000000 | |
| Stellar_Type(1)>MT | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | ... | 8.000000 | 8.000000 | 8.000000 | 13.000000 | 13.000000 | 7.000000 | 13.000000 | 13.000000 | 7.000000 | 13.000000 | |
| Stellar_Type(2)<MT | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | ... | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 2.000000 | 8.000000 | 1.000000 | 2.000000 | |
| Stellar_Type(2)>MT | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | ... | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 7.000000 | 13.000000 | 1.000000 | 7.000000 | |
| Time<MT | Myr | 18.796846 | 18.800936 | 18.804943 | 18.808871 | 18.812720 | 18.816491 | 18.820188 | 18.823810 | 18.827360 | ... | 18.885181 | 18.887504 | 18.889780 | 18.892010 | 19.892010 | 7.252029 | 12.304917 | 13.459892 | 6.795883 | 11.419296 |
| Time>MT | Myr | 18.796846 | 18.800936 | 18.804943 | 18.808871 | 18.812720 | 18.816491 | 18.820188 | 18.823810 | 18.827360 | ... | 18.885181 | 18.887504 | 18.889780 | 18.892010 | 19.892010 | 7.252029 | 12.304917 | 13.459892 | 6.795883 | 11.419296 |
| Zeta_Hurley(1) | - | 0.437277 | 0.345061 | 0.345982 | 0.347154 | 0.348360 | 0.349597 | 0.350865 | 0.352166 | 0.353500 | ... | 0.388655 | 0.390922 | 0.393248 | 0.000000 | 0.000000 | 0.300000 | 0.000000 | 0.000000 | 0.300000 | 0.000000 |
| Zeta_Hurley(2) | - | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | ... | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.000000 | 0.300000 | 0.300000 |
| Zeta_Hurley_He(1) | - | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| Zeta_Hurley_He(2) | - | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | ... | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.000000 | 0.000000 | 0.300000 | 0.000000 |
| Zeta_SoberMan_He(1) | - | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| Zeta_SoberMan_He(2) | - | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | ... | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | 0.000000 | 0.000000 | -0.333333 | 0.000000 |
| Zeta_Soberman(1) | - | 0.859915 | 0.483761 | 0.488665 | 0.494846 | 0.501124 | 0.507484 | 0.513929 | 0.520460 | 0.527078 | ... | 0.681138 | 0.690121 | 0.699247 | 0.000000 | 0.000000 | -0.333333 | 0.000000 | 0.000000 | -0.333333 | 0.000000 |
| Zeta_Soberman(2) | - | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | ... | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | 0.000000 | -0.333333 | -0.333333 |
32 rows × 39 columns
getEventHistory¶
Often, it is useful to quickly retrieve an overview of the event history of a binary, including all mass transfer, common envelope, and supernova events. For this, we use getEventHistory
[10]:
seeds, events = getEventHistory(Data)
for ii, seed in enumerate(seeds):
print(seed, events[ii])
1636090318 []
1636090389 []
1636091116 []
101048 [('MT', 7.252028833917178, 2, 1, True, False, False), ('SN', 8.223359358765522, 8, 13, 1, False), ('MT', 12.304916903390055, 13, 2, False, True, True), ('MT', 13.459891998705004, 13, 8, False, True, False), ('SN', 13.459891998705004, 8, 13, 2, False)]
102660 [('MT', 6.795882556458199, 2, 1, True, False, False), ('SN', 7.702682726277025, 8, 13, 1, False), ('MT', 11.419295515325874, 13, 2, False, True, True), ('SN', 12.561084237007242, 8, 13, 2, False)]
getEventHistory takes the h5file as input, and returns an array of the seeds processed as well as the major events for that seed. The format for events depends on the event type. Currently, we only include supernova and mass transfer events, but mass transfer events include a flag for whether the system underwent CEE.
For MT events, it is (‘MT’, time, stellarType1, stellarType2, isRlof1, isRlof2, isCEE)
For SN events, it is (‘SN’, time, stellarTypeProgenitor, stellarTypeRemnant, whichIsProgenitor, isUnbound)
There is also an optional argument exclude_null which defaults to False. If True, it will skip systems which undergo no events of interest (which may speed up large runs).
A useful function that builds off of getEventHistory is getEventStrings, which collects the event information into a succint string, which may be easier to read (once you get used to the syntax).
The syntax for the event strings takes the following convention:
For MT events:
P>S, P<S, or P=S
where P is primary type, S is secondary type, and >, < is RLOF (1->2 or 1<-2) or = for CEE
For SN events:
P*SR for star1 the SN progenitor, or
R*SP for star2 the SN progenitor,
where P is progenitor type, R is remnant type, S is state (I for intact, U for unbound)
Event strings for the same seed are ordered chronologically and separated by the undesrcore character _
[11]:
eventStrings = getEventStrings(allEvents=events)
for ii in range(5):
print(seeds[ii], eventStrings[ii])
1636090318 NA
1636090389 NA
1636091116 NA
101048 2>1_8*I13_13=2_13<8_13*I8
102660 2>1_8*I13_13=2_13*I8
[ ]:
2. Slicing the data¶
Since the random seed is unique and constant for a given binary, the properties and events of the binary system can be recovered by looking at its seed across different output categories.
Here we introduce the basics of manipulating the data using the seeds. We provide an example on how we get the initial parameters of systems that ended up forming double compact objects.
Naively, we might try to use For Loops with Conditions to extract systems of interest to a list. However, this can potentially be computationally expensive.
Here we present a method to more efficiently ‘slice’ the data using numpy and boolean masks. These are slightly more involved but are computationally quick and use intuitive logic.
Question: What were the initial total masses of the double compact objects?¶
[12]:
def calculateTotalMassesNaive(pathData=None):
Data = h5.File(pathToData)
totalMasses = []
# Retrive the categories
SPs = Data['BSE_System_Parameters']
DCs = Data['BSE_Double_Compact_Objects']
# For syntax see section 1
# Extract parameters of interest
seedsDC = DCs['SEED'][()]
seedsSP = SPs['SEED'][()]
m1Zams = SPs['Mass@ZAMS(1)'][()]
m2Zams = SPs['Mass@ZAMS(2)'][()]
for dcSeed in seedsDC:
for seedIndex in range(len(seedsSP)):
spSeed = seedsSP[seedIndex]
if spSeed == dcSeed:
m1 = m1Zams[seedIndex]
m2 = m2Zams[seedIndex]
mTot = m1 + m2
totalMasses.append(mTot)
Data.close()
return totalMasses
[13]:
# calculate function run time
start = time.time()
mTotOld = calculateTotalMassesNaive(pathData=pathToData)
end = time.time()
timeDiffNaive = end-start
print('%s seconds, using for loops.' %(timeDiffNaive))
0.0061054229736328125 seconds, using for loops.
I) Optimizing the above loop¶
a - Use built-in numpy routines¶
Numpy arrays can make use of a powerful library of optimization tools which allow the user to bypass computationally heavy for loops.
For example, we can speed up the calculation of the element-wise sum of two arrays with:
[14]:
SPs = Data['BSE_System_Parameters']
m1Zams = SPs['Mass@ZAMS(1)'][()]
m2Zams = SPs['Mass@ZAMS(2)'][()]
mTotalAllSystems = np.add(m1Zams, m2Zams)
b - Use boolean masks in a single file¶
Where previously we put the condition in an if statement nested within a for loop, now we again make use of boolean masks to filter out the undesired elements.
The boolean array must have the same length as the input array.
[15]:
# Create a boolean array from the total mass array which is True
# if the total mass of the corrresponding system is less than 40.
maskMTotLessThan40 = (mTotalAllSystems <= 40)
Crucially, you can apply this mask to all other columns in the same file because, by construction, they all have the same length.
[16]:
# seeds of systems with total mass below 40
seedsMtotBelow40 = seedsSP[maskMTotLessThan40]
Note that this works because the order of the two columns (seeds and total masses) are the same.
For example, the total mass of the third system entry corresponds to the seed at the third system entry.
II) Use seeds as masks between files¶
Example 1¶
Before we continue it is useful to understand how the COMPAS-popsynth printing works.
Each simulated system will be initialized only once and so will have only one line in the BSE_System_Parameters file. However, lines in BSE_RLOF are created whenever a system goes through a mass transfer event, which might happen multiple times for a single system, or potentially not at all. Similarly, in the BSE_Supernovae file, you will find at most two lines per system, but possibly none. BSE_Double_Compact_Objects lines are printed only when the final system is intact and
composed of either Neutron Stars or Black Holes, which is a rare event that happens at most once per system.
For this reason, it is in general not the case that the system on line \(n\) of one file corresponds will match the system on line \(n\) of another file.
In order to match systems across files, we need to extract the seeds of desired systems from one file, and apply them as a mask in the other file.
[17]:
# Example: calculate the primary ZAMS mass of systems which become DCOs (Double Compact Objects)
seedsSP = SPs['SEED'][()]
seedsDC = DCs['SEED'][()]
m1Zams = SPs['Mass@ZAMS(1)'][()]
# Calculate mask for which elements of seedsSP are found in seedsDC
# - see numpy.in1d documentation for details
mask = np.in1d(seedsSP, seedsDC)
print(mask)
print(seedsDC)
print(seedsSP[mask])
print(m1Zams[mask])
print("The occurence rate of DCOs is {}/{}".format(sum(mask), len(mask)))
[False True True True True]
[1636090389 1636091116 101048 102660]
[1636090389 1636091116 101048 102660]
[ 76.555207 113.445902 24.94177 26.74675 ]
The occurence rate of DCOs is 4/5
[18]:
printCompasDetails(DCs, [1636090389, 1636091089, 1636091116])
[18]:
| SEED | (units) | 1636090389 | 1636091116 |
|---|---|---|---|
| Coalescence_Time | Myr | 447.467040 | 5.688494e+07 |
| Eccentricity@DCO | - | 0.346775 | 1.280522e-01 |
| Mass(1) | Msol | 9.883289 | 1.332088e+01 |
| Mass(2) | Msol | 2.315392 | 1.332038e+01 |
| Merges_Hubble_Time | State | 1.000000 | 0.000000e+00 |
| Recycled_NS(1) | Event | 0.000000 | 0.000000e+00 |
| Recycled_NS(2) | Event | 0.000000 | 0.000000e+00 |
| SemiMajorAxis@DCO | AU | 0.027903 | 9.694038e-01 |
| Stellar_Type(1) | - | 14.000000 | 1.400000e+01 |
| Stellar_Type(2) | - | 13.000000 | 1.400000e+01 |
| Time | Myr | 6.750712 | 5.321871e+00 |
Optimized loop¶
[19]:
def calculateTotalMassesOptimized(pathData=None):
Data = h5.File(pathToData)
totalMasses = []
# Retrive the categories
SPs = Data['BSE_System_Parameters']
DCs = Data['BSE_Double_Compact_Objects']
# For syntax see section 1
# Extract parameters of interest
seedsDC = DCs['SEED'][()]
seedsSP = SPs['SEED'][()]
m1Zams = SPs['Mass@ZAMS(1)'][()]
m2Zams = SPs['Mass@ZAMS(2)'][()]
mZamsTot = np.add(m1Zams, m2Zams)
maskSeedsBecameDCO = np.in1d(seedsSP, seedsDC)
mZamsTotOfDCOs = mZamsTot[maskSeedsBecameDCO]
Data.close()
return mZamsTotOfDCOs
[20]:
# calculate function run time
start = time.time()
mTotNew = calculateTotalMassesOptimized(pathData=pathToData)
end = time.time()
timeDiffOptimized = end-start
# calculate number of Double Compact Objects
nrDCOs = len(seedsDC)
print('Compare')
print('%s seconds, using For Loops.' %(timeDiffNaive))
print('%s seconds, using Optimizations.' %(timeDiffOptimized))
print('Using %s DCO systems' %(nrDCOs))
Compare
0.0061054229736328125 seconds, using For Loops.
0.003986835479736328 seconds, using Optimizations.
Using 4 DCO systems
Note: The time difference will depend heavily on the number of systems under investigation, as well as the number of bypassed For Loops. If you used the path to the pre-generated tutorial data set (with few, intentionally specified systems), you should see very little improvement.
[21]:
# Test that the two arrays are in fact identical
print(np.array_equal(mTotOld, mTotNew))
True
Note: the above loop can easily be expanded with more conditions.
If you do not want all the DCO initial total masses but only of the double neutron stars, then you just need to apply another mask to the seedsDC.
[22]:
def calculateTotalMassesDNS(pathToData=None):
Data = h5.File(pathToData)
totalMasses = []
SPs = Data['BSE_System_Parameters']
DCs = Data['BSE_Double_Compact_Objects']
seedsDC = DCs['SEED'][()]
stype1 = DCs['Stellar_Type(1)'][()]
stype2 = DCs['Stellar_Type(2)'][()]
dcMaskDNS = (stype1 == 13) & (stype2 == 13)
seedsDNS = seedsDC[dcMaskDNS]
# Get info from ZAMS
seedsSP = SPs['SEED'][()]
m1Zams = SPs['Mass@ZAMS(1)'][()]
m2Zams = SPs['Mass@ZAMS(2)'][()]
mZamsTot = np.add(m1Zams, m2Zams)
spMaskDNS = np.in1d(seedsSP, seedsDNS)
mZamsTotDNS = mZamsTot[spMaskDNS]
Data.close()
return mZamsTotDNS
[23]:
# calculate function run time
start = time.time()
mTotDNS = calculateTotalMassesDNS(pathToData=pathToData)
end = time.time()
timeDiffDNS = end-start
# calculate number of DNS systems
nrDNSs = len(mTotDNS)
print('%s seconds for all %s DNS systems.' %(timeDiffDNS, nrDNSs))
0.004500150680541992 seconds for all 2 DNS systems.
The printCompasDetails function can also optionally take a mask as argument. This is especially useful for those output categories which have multiple events for a single seed. Using both seeds and mask inputs can help to extract a specifc type of event from several for the given seeds.
The mask array must have the same length as the data arrays for the given category.
[24]:
# Example: Want to investigate CEE events for seed 1636090318
printCompasDetails(MTs, 1636090318)
[24]:
| SEED | (units) | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | ... | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 | 1636090318 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| CEE>MT | State | 1.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 1.000000 |
| Eccentricity<MT | - | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.743907 |
| Eccentricity>MT | - | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.743907 | 0.000000 |
| MT_Event_Counter | Count | 1.000000 | 2.000000 | 3.000000 | 4.000000 | 5.000000 | 6.000000 | 7.000000 | 8.000000 | 9.000000 | ... | 24.000000 | 25.000000 | 26.000000 | 27.000000 | 28.000000 | 29.000000 | 30.000000 | 31.000000 | 32.000000 | 33.000000 |
| Mass(1)<MT | Msol | 11.773666 | 3.943185 | 4.952215 | 4.957205 | 4.957273 | 4.957274 | 4.957274 | 4.957275 | 4.957275 | ... | 4.957278 | 4.957278 | 4.957279 | 4.957279 | 4.957279 | 4.957279 | 4.957279 | 4.957279 | 4.957279 | 1.360533 |
| Mass(1)>MT | Msol | 3.943185 | 4.952215 | 4.957205 | 4.957273 | 4.957274 | 4.957274 | 4.957275 | 4.957275 | 4.957275 | ... | 4.957278 | 4.957279 | 4.957279 | 4.957279 | 4.957279 | 4.957279 | 4.957279 | 4.957279 | 1.360533 | 1.360533 |
| Mass(2)<MT | Msol | 2.403020 | 2.403020 | 1.393990 | 1.389000 | 1.388932 | 1.388931 | 1.388931 | 1.388930 | 1.388930 | ... | 1.388927 | 1.388927 | 1.388926 | 1.388926 | 1.388926 | 1.388926 | 1.388926 | 1.388926 | 1.388925 | 1.388925 |
| Mass(2)>MT | Msol | 2.403020 | 1.393990 | 1.389000 | 1.388932 | 1.388931 | 1.388931 | 1.388930 | 1.388930 | 1.388930 | ... | 1.388927 | 1.388926 | 1.388926 | 1.388926 | 1.388926 | 1.388926 | 1.388926 | 1.388925 | 1.388925 | 1.388925 |
| RLOF(1)<MT | Event | 0.000000 | 1.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| RLOF(1)>MT | State | 1.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| RLOF(2)<MT | Event | 0.000000 | 0.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | ... | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 |
| RLOF(2)>MT | State | 0.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | ... | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 |
| Radius(1)<MT | Rsol | 569.606809 | 0.614182 | 0.500045 | 0.505819 | 0.511921 | 0.518248 | 0.524816 | 0.531644 | 0.538753 | ... | 0.702361 | 0.719599 | 0.738139 | 0.758126 | 0.779723 | 0.803117 | 0.828516 | 0.856161 | 0.886326 | 0.000014 |
| Radius(1)>MT | Rsol | 0.597875 | 0.494328 | 0.499924 | 0.505817 | 0.511921 | 0.518248 | 0.524816 | 0.531644 | 0.538753 | ... | 0.702361 | 0.719599 | 0.738139 | 0.758126 | 0.779723 | 0.803116 | 0.828516 | 0.856161 | 0.000014 | 0.000014 |
| Radius(2)<MT | Rsol | 1.698892 | 1.698897 | 1.357432 | 1.353888 | 1.353840 | 1.353839 | 1.353840 | 1.353840 | 1.353840 | ... | 1.353843 | 1.353844 | 1.353844 | 1.353844 | 1.353844 | 1.353844 | 1.353845 | 1.353845 | 1.353845 | 1.353968 |
| Radius(2)>MT | Rsol | 1.698892 | 1.357432 | 1.353887 | 1.353839 | 1.353839 | 1.353839 | 1.353839 | 1.353840 | 1.353840 | ... | 1.353843 | 1.353843 | 1.353844 | 1.353844 | 1.353844 | 1.353844 | 1.353844 | 1.353845 | 1.353845 | 1.353968 |
| SemiMajorAxis<MT | Rsol | 1079.189435 | 2.594059 | 4.880604 | 4.905801 | 4.906144 | 4.906150 | 4.906152 | 4.906153 | 4.906154 | ... | 4.906171 | 4.906172 | 4.906173 | 4.906173 | 4.906174 | 4.906175 | 4.906176 | 4.906177 | 4.906178 | 2.936388 |
| SemiMajorAxis>MT | Rsol | 2.594059 | 4.880604 | 4.905801 | 4.906144 | 4.906150 | 4.906152 | 4.906153 | 4.906154 | 4.906156 | ... | 4.906172 | 4.906173 | 4.906173 | 4.906174 | 4.906175 | 4.906176 | 4.906177 | 4.906178 | 2.936388 | 0.751990 |
| Stellar_Type(1)<MT | 5.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | ... | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 13.000000 | |
| Stellar_Type(1)>MT | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | ... | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 8.000000 | 13.000000 | 13.000000 | |
| Stellar_Type(2)<MT | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | ... | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | |
| Stellar_Type(2)>MT | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | ... | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | |
| Time<MT | Myr | 18.796846 | 18.800936 | 18.804943 | 18.808871 | 18.812720 | 18.816491 | 18.820188 | 18.823810 | 18.827360 | ... | 18.872837 | 18.875407 | 18.877925 | 18.880392 | 18.882811 | 18.885181 | 18.887504 | 18.889780 | 18.892010 | 19.892010 |
| Time>MT | Myr | 18.796846 | 18.800936 | 18.804943 | 18.808871 | 18.812720 | 18.816491 | 18.820188 | 18.823810 | 18.827360 | ... | 18.872837 | 18.875407 | 18.877925 | 18.880392 | 18.882811 | 18.885181 | 18.887504 | 18.889780 | 18.892010 | 19.892010 |
| Zeta_Hurley(1) | - | 0.437277 | 0.345061 | 0.345982 | 0.347154 | 0.348360 | 0.349597 | 0.350865 | 0.352166 | 0.353500 | ... | 0.378139 | 0.380137 | 0.382187 | 0.384289 | 0.386444 | 0.388655 | 0.390922 | 0.393248 | 0.000000 | 0.000000 |
| Zeta_Hurley(2) | - | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | ... | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 |
| Zeta_Hurley_He(1) | - | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| Zeta_Hurley_He(2) | - | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | ... | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 | 0.300000 |
| Zeta_SoberMan_He(1) | - | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| Zeta_SoberMan_He(2) | - | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | ... | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 |
| Zeta_Soberman(1) | - | 0.859915 | 0.483761 | 0.488665 | 0.494846 | 0.501124 | 0.507484 | 0.513929 | 0.520460 | 0.527078 | ... | 0.638267 | 0.646580 | 0.655020 | 0.663591 | 0.672296 | 0.681138 | 0.690121 | 0.699247 | 0.000000 | 0.000000 |
| Zeta_Soberman(2) | - | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | ... | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 | -0.333333 |
32 rows × 34 columns
There is too much information here to be useful in this case, so we apply a mask to filter out events we’re not interested in
[25]:
maskCEE = MTs['CEE>MT'][()] == 1 # systems which undergo Common Envelope Evolution
printCompasDetails(MTs, 1636090318, mask=maskCEE)
[25]:
| SEED | (units) | 1636090318 | 1636090318 |
|---|---|---|---|
| CEE>MT | State | 1.000000 | 1.000000 |
| Eccentricity<MT | - | 0.000000 | 0.743907 |
| Eccentricity>MT | - | 0.000000 | 0.000000 |
| MT_Event_Counter | Count | 1.000000 | 33.000000 |
| Mass(1)<MT | Msol | 11.773666 | 1.360533 |
| Mass(1)>MT | Msol | 3.943185 | 1.360533 |
| Mass(2)<MT | Msol | 2.403020 | 1.388925 |
| Mass(2)>MT | Msol | 2.403020 | 1.388925 |
| RLOF(1)<MT | Event | 0.000000 | 0.000000 |
| RLOF(1)>MT | State | 1.000000 | 0.000000 |
| RLOF(2)<MT | Event | 0.000000 | 1.000000 |
| RLOF(2)>MT | State | 0.000000 | 1.000000 |
| Radius(1)<MT | Rsol | 569.606809 | 0.000014 |
| Radius(1)>MT | Rsol | 0.597875 | 0.000014 |
| Radius(2)<MT | Rsol | 1.698892 | 1.353968 |
| Radius(2)>MT | Rsol | 1.698892 | 1.353968 |
| SemiMajorAxis<MT | Rsol | 1079.189435 | 2.936388 |
| SemiMajorAxis>MT | Rsol | 2.594059 | 0.751990 |
| Stellar_Type(1)<MT | 5.000000 | 13.000000 | |
| Stellar_Type(1)>MT | 8.000000 | 13.000000 | |
| Stellar_Type(2)<MT | 1.000000 | 1.000000 | |
| Stellar_Type(2)>MT | 1.000000 | 1.000000 | |
| Time<MT | Myr | 18.796846 | 19.892010 |
| Time>MT | Myr | 18.796846 | 19.892010 |
| Zeta_Hurley(1) | - | 0.437277 | 0.000000 |
| Zeta_Hurley(2) | - | 0.300000 | 0.300000 |
| Zeta_Hurley_He(1) | - | 0.000000 | 0.000000 |
| Zeta_Hurley_He(2) | - | 0.300000 | 0.300000 |
| Zeta_SoberMan_He(1) | - | 0.000000 | 0.000000 |
| Zeta_SoberMan_He(2) | - | -0.333333 | -0.333333 |
| Zeta_Soberman(1) | - | 0.859915 | 0.000000 |
| Zeta_Soberman(2) | - | -0.333333 | -0.333333 |
[ ]:
Example 2¶
The previous example uses the fact that both SystemParameters and DoubleCompactObjects only print at most one line per system. However, as mentioned above, events such as supernovae or common envelopes might happen multiple times to a given system, and as a result there would be multiple occurences of a given seed in the relevant file.
To account for this, we will need to modify the previous method. Consider again the 4 seeds of the previous example. Both 2 and 4 formed a DCO and hence both stars in these binaries went SN. Seeds 1 and 3 are low mass stars hence they did not go SN. (Note that we do not specify the companion masses for any of these systems, but for simplicity we assume that the companions to 1 and 3 are also sufficiently low mass to not produce a supernova). The SN file prints one line per SN and therefore seeds 2 and 4 appear twice each.
Imagine you want the primary masses of systems that experienced at any point a core collapse supernova (CCSN). We’ll reuse our mock data, with additional information about the types of SN which occured in each star. Here, PPISN refers to Pulsational Pair Instability Supernovae.
[26]:
printCompasDetails(SNe)
[26]:
| SEED | (units) | 1636090318 | 1636090389 | 1636090389 | 1636091116 | 1636091116 | 101048 | 101048 | 102660 | 102660 |
|---|---|---|---|---|---|---|---|---|---|---|
| Applied_Kick_Magnitude(SN) | kms^-1 | 505.642361 | 45.113582 | 50.969609 | 16.648911 | 23.878749 | 76.914889 | 66.973808 | 184.916045 | 180.417034 |
| ComponentSpeed(CP) | kms^-1 | 386.014275 | 5.491568 | 128.586076 | 2.041058 | 14.649187 | 16.926151 | 177.401632 | 28.402032 | 129.864312 |
| ComponentSpeed(SN) | kms^-1 | 386.014275 | 5.491568 | 128.586076 | 2.041058 | 14.649187 | 16.926151 | 177.401632 | 28.402032 | 129.864312 |
| Drawn_Kick_Magnitude(SN) | kms^-1 | 549.417344 | 217.320105 | 50.969609 | 435.503588 | 624.150245 | 100.550626 | 66.973808 | 249.009476 | 225.718097 |
| Eccentricity | - | 0.743907 | 0.4113 | 0.346775 | 0.050379 | 0.128052 | 0.707907 | 0.295882 | 0.859918 | 0.474411 |
| Eccentricity<SN | - | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| Experienced_RLOF(SN) | Event | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
| Fallback_Fraction(SN) | - | 0.079675 | 0.79241 | 0.0 | 0.961771 | 0.961742 | 0.235063 | 0.0 | 0.257394 | 0.200698 |
| Is_Hydrogen_Poor(SN) | State | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
| Kick_Magnitude(uK) | - | 1.018117 | 0.254925 | 0.059579 | 0.104508 | 0.160635 | 0.594548 | 0.083929 | 1.736729 | 0.463334 |
| MT_Donor_Hist(SN) | - | b'5 ' | b'2 ' | b'2-7-8 ' | b'2 ' | b'2 ' | b'2 ' | b'2-8 ' | b'2 ' | b'2 ' |
| Mass(CP) | Msol | 1.388925 | 60.96837 | 9.883289 | 96.197042 | 13.320884 | 20.646428 | 2.008356 | 21.113558 | 2.155738 |
| Mass(SN) | Msol | 1.360533 | 9.879445 | 2.315392 | 13.320862 | 13.320378 | 2.008319 | 1.587026 | 2.155738 | 1.77577 |
| Mass_CO_Core@CO(SN) | Msol | 2.698725 | 9.331406 | 5.179326 | 10.695529 | 10.695299 | 4.131764 | 3.25185 | 4.435399 | 3.515565 |
| Mass_Core@CO(SN) | Msol | 2.698725 | 9.331406 | 5.179326 | 10.695529 | 10.695299 | 4.131764 | 3.25185 | 4.435399 | 3.515565 |
| Mass_He_Core@CO(SN) | Msol | 4.957279 | 12.227034 | 5.179326 | 13.898675 | 13.898547 | 5.600205 | 3.25185 | 5.978744 | 4.848981 |
| Mass_Total@CO(SN) | Msol | 4.957279 | 12.227034 | 5.179326 | 13.898675 | 13.898547 | 5.600205 | 3.25185 | 5.978744 | 4.848981 |
| Orb_Velocity<SN | kms^-1 | 496.644721 | 176.968272 | 855.489186 | 159.307182 | 148.65178 | 129.367001 | 797.979295 | 106.473722 | 389.389086 |
| SN_Kick_Mean_Anomaly(SN) | - | 5.189746 | 1.824014 | 5.662449 | 0.940618 | 1.977927 | 5.294177 | 2.189489 | 3.100337 | 5.305596 |
| SN_Kick_Phi(SN) | - | 0.648572 | -0.663334 | 0.390192 | -0.17818 | 0.598305 | -0.748282 | 0.095727 | 1.205367 | -0.340507 |
| SN_Kick_Theta(SN) | - | 3.478068 | 3.615591 | 0.542629 | 1.748962 | 0.997198 | 5.110205 | 5.717921 | 5.544482 | 3.657304 |
| SN_Type(SN) | - | 1 | 1 | 16 | 1 | 1 | 1 | 16 | 1 | 1 |
| SemiMajorAxis | Rsol | 2.936388 | 725.244673 | 6.000134 | 870.531899 | 208.452875 | 981.529584 | 2.230194 | 2880.792791 | 12.763265 |
| SemiMajorAxis<SN | Rsol | 4.906178 | 445.670493 | 3.924571 | 827.219197 | 234.88739 | 299.052528 | 1.575217 | 455.703044 | 8.809354 |
| Stellar_Type(CP) | - | 1 | 1 | 14 | 1 | 14 | 1 | 13 | 1 | 13 |
| Stellar_Type(SN) | - | 13 | 14 | 13 | 14 | 14 | 13 | 13 | 13 | 13 |
| Stellar_Type_Prev(SN) | - | 8 | 8 | 8 | 8 | 8 | 8 | 8 | 8 | 8 |
| Supernova_State | State | 1 | 1 | 2 | 1 | 2 | 1 | 2 | 1 | 2 |
| SystemicSpeed | kms^-1 | 386.014275 | 5.491568 | 128.586076 | 2.041058 | 14.649187 | 16.926151 | 177.401632 | 28.402032 | 129.864312 |
| Time | Myr | 18.89201 | 4.252851 | 6.750712 | 3.872097 | 5.321871 | 8.223359 | 13.459892 | 7.702683 | 12.561084 |
| Unbound | State | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
[27]:
# Example: get the primary ZAMS masses of systems which experience 2 CCSNe before becoming a DCO
seedsSP = SPs['SEED'][()]
seedsSN = SNe['SEED'][()]
seedsDC = DCs['SEED'][()]
snType = SNe['SN_Type(SN)'][()]
m1Zams = SPs['Mass@ZAMS(1)'][()]
# Note: the SN_Type(SN) parameter maps integers to SN types:
snTypeDict = {
1: 'CCSN',
2: 'ECSN',
16: 'USSN'
} # this dictionary is illustrative, but not explicitly used here
# Determine which seeds experienced at least 1 CCSN
maskCCSN = snType == 1
seedsCCSN, countsCCSN = np.unique(seedsSN[maskCCSN], return_counts=True)
# Seeds with 2 CCSNe will have a countsCCSN value of 2
seedsDoubleCCSN = seedsCCSN[countsCCSN == 2]
# Make a mask for SPs using the seeds aquired above
maskDoubleCcsnSP = np.in1d(seedsSP, seedsDoubleCCSN)
m1ZamsDoubleCcsn = m1Zams[maskDoubleCcsnSP]
print("Primary ZAMS masses for systems which undergo 2 CCSNe before becoming a DCO are: ",m1ZamsDoubleCcsn)
Primary ZAMS masses for systems which undergo 2 CCSNe before becoming a DCO are: [113.445902 26.74675 ]
[28]:
# Always remember to close your data file
Data.close()
3. Visualizing the data¶
Although math is the fundamental basis of physics and astrophysics, we cannot always easily convert numbers and equations into a coherent picture. Plotting is therefore a vital tool in bridging the gap between raw data and a deeper scientific understanding.
Disclaimer:
There are many ways to make the same plot in matplotlib and there are many ways to bin your data. Often, there is no “best” way to display data in a plot, and the message conveyed can be heavily dependent on the context of the data as well as asthetic plotting decisions.
For example, in histograms, as we discuss below, the relatively subjective choice of bin size can significantly affect the interpretation of the results. It is important to be aware of when and how we make these choices and to try to reduce any unintended bias.
** Example: inspect the component masses of Double Compact Objects**
In the example below, we use the following conventions:
1 - We deliberately choose to use the matplotlib.pyplot.subplots routine even when creating a single figure (as opposed to using pyplot.plot). This is because many online forums (e.g Stackoverflow) use this syntax. Furthermore, this means you do not have to learn two different types of syntax when creating either a single or multiple panel figure.
2 - We choose to do the binning within the numpy/array environment instead of with inbuilt functions such as plt.hist / axes.hist. The reason is that you have more control over what you do, such as custom normalization (using rates, weights, pdf, etc.). It also forces you to have a deeper understanding of what you are calculating, and allows you to check intermediate steps with print statements. Once you know how to bin your data this way you can also easily expand these routines for more complicated plots (2D binning).
Note: for this exercise, we recommend running your own simulation of at least 100,000 binaries in order to have a sufficient number of DCOs to have an interesting plot. We use the default tutorial data here for illustrative purposes, and because such a large data file would be too large to store on github.
Get some data to plot¶
[30]:
pathToData = 'COMPAS_Tutorial_Output.h5'
Data = h5.File(pathToData)
print(list(Data.keys()))
DCs = Data['BSE_Double_Compact_Objects']
m1 = DCs['Mass(1)'][()]
m2 = DCs['Mass(2)'][()]
mTot = np.add(m1, m2)
Data.close()
['BSE_Common_Envelopes', 'BSE_Double_Compact_Objects', 'BSE_RLOF', 'BSE_Supernovae', 'BSE_System_Parameters', 'Run_Details']
Plot histogram and CDF of data on left, and component mass scatter plot on the right¶
[54]:
fig, axes = plt.subplots(ncols=2, figsize=(15,8))
largefontsize = 30
smallfontsize = 20
# Histogram
ax1 = axes[0]
bins = np.linspace(0, max(mTot), 21) # use 20 bins, up to the maximum total DCO mass
ax1.hist(mTot, bins=bins)
ax1.set_ylabel('Histogram counts', fontsize=smallfontsize)
ax1.set_xlabel('Total mass', fontsize=smallfontsize)
ax1.set_title('Total mass at DCO formation', fontsize=largefontsize)
# CDF
ax2 = ax1.twinx()
cdf_xvalues = np.cumsum(np.sort(mTot))
np.insert(cdf_xvalues, 0, 0) # insert a 0 at the front of the array
cdf_yvalues = np.linspace(0, ax2.get_ylim()[1], len(cdf_xvalues))
ax2.plot(cdf_xvalues, cdf_yvalues)
ax2.set_ylabel('CDF values', rotation=270, labelpad=15)
ax2.grid(False)
# Scatter plot
ax3 = axes[1]
ax3.scatter(m1, m2)
ax3.set_xlabel(r'M1 [$M_\odot$]', fontsize=smallfontsize)
ax3.set_ylabel(r'M2 [$M_\odot$]', fontsize=smallfontsize)
ax3.set_title('Component Masses', fontsize=largefontsize)
fig.tight_layout()
