<jupyter_start><jupyter_text>Demo: Chirp mass distribution of binary black holes (BBHs)This notebook reproduces Figure 10 of the COMPAS Paper (_arXiv_: ; Riley et al., in prep.) The data can be reproduced by running the attached `stroopwafel_interface.py`, that reads in the `pythonSubmit.py`This file contains the output of 10,000,000 binaries evolved using COMPAS v02.19.04, using adaptive importance sampling , sampling from a metallicity uniform in $\log(Z) \in [10^{-4},0.03]$. More details can be found in `Run_Details.txt`.The rates were calculated by running ```FastCosmicIntegration.py```  from COMPAS's post-processing tools, with the following flags altered from their default values:```:::bash    python FastCosmicIntegration.py  --path /n/de_mink_lab/Users/lvanson/CompasOutput/v02.19.04/N1e7_Fid/fWR1.0coolWind1.0//output/ --filename COMPAS_Output_wWeights.h5 --mu0 0.035 --muz -0.23 --sigma0 0.39 --sigmaz 0.0 --alpha 0.0 --weight mixture_weight --zstep 0.01 --sens O1 --m1min 10. --aSF 0.01 --bSF 2.77 --cSF 2.9 --dSF 4.7 > $OUT_DIR/masterfolder/CosmicIntegration/COMPAS_CI1.log ```For a more detailed explanation of the workings of COMPAS's post-processing tools, please work through the Jupyter Notebook tutorials in the `postProcessing` directory.  Paths and libraries* If you would like to use this script for your own COMPAS Output, you can change `fileLoc`.<jupyter_code>import numpy as np
import h5py as h5
from astropy.table import Table, Column

import astropy.units as u
from astropy import constants as const

import matplotlib.pyplot as plt
from matplotlib import ticker, cm

from scipy import stats


######################################
fileLoc    = '/Volumes/TOSHIBA_EXT/CompasOutput/v02.19.04/N1e7_Fid_forCOMPASmethod/output/'

######################################
##PLOT setttings
plt.rc('font', family='serif')

from matplotlib import rc
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'stix'
matplotlib.rcParams['font.family'] = 'STIXGeneral'
rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
rc('text', usetex=True)
fsize, SMALL_SIZE, MEDIUM_SIZE, BIGGER_SIZE = 30,25,25,30
for obj in ['axes','xtick','ytick']:
    plt.rc(obj, labelsize=MEDIUM_SIZE)          # controls default text sizes
for obj in ['figure','axes']:
    plt.rc(obj, titlesize=BIGGER_SIZE)    # fontsize of the tick labels
plt.rc('font', size=MEDIUM_SIZE)          # controls default text sizes
plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize<jupyter_output><empty_output><jupyter_text>Read the data<jupyter_code>################################################
## Open hdf5 file
File        = h5.File(fileName+'COMPAS_Output_wWeights.h5' ,'r')

# The group that contains the rate data
rate_key = 'Rates_mu00.035_muz-0.23_alpha0.0_sigma00.39_sigmaz0.0'
################################################
## Read merger rate related data
DCO_mask                  = File[rate_key]['DCOmask'][()] # Mask from DCO to merging BBH 
redshifts                 = File[rate_key]['redshifts'][()] # Redshifts at which the rates were calculated
O1_detection_rate         = File[rate_key]['detection_rateO1'][()] # detection rate for O1 sensitivity 
# intrinsic_rate_density    = File[rate_key]['merger_rate'][()]  # Intrinsic BBH merger rate density


# Read the DCO masses, Select only merging BBHs
Mass1              = File['DoubleCompactObjects']['Mass(1)'][DCO_mask]
Mass2              = File['DoubleCompactObjects']['Mass(2)'][DCO_mask]
M_moreMassive      = np.maximum(Mass1, Mass2)
M_lessMassive      = np.minimum(Mass1, Mass2)


File.close()<jupyter_output><empty_output><jupyter_text>Set values<jupyter_code>######################################
def Mchirp(m1, m2):
    chirp_mass = np.divide(np.power(np.multiply(m1, m2), 3./5.), np.power(np.add(m1, m2), 1./5.))
    return chirp_mass    
# Calculate the chrip mass
Mchirp = Mchirp(M_moreMassive, M_lessMassive)

# Use the detection rates at O1 sensitivity
ratePerSys_observed = np.sum(O1_detection_rate, axis=1)
mask = Mchirp < 100 # Set maximum plotted chirp mass to be 100 Msun

# Set grid in chirp mass
binsMchirp  = np.linspace(0, 40, 20)
binwidth    = np.diff(binsMchirp)
center_bins = (binsMchirp[1:] + binsMchirp[:-1]) / 2.0

# yvalues, _ = np.histogram(Mchirp[mask], bins=binsMchirp, weights=ratePerSys_observed[mask])
# dydMchirp = np.divide(yvalues, dMchirp)

#################################
# The data that you will plot
x_vals  = Mchirp
Weights = ratePerSys_observed
Bins    = binsMchirp

#################################
# Would you like to bootstrap your data?
# Warning, this is slow
bootstrap  = True
bootstraps =  50

#################################
# What colour should your plot be
colors  = 'grey'
#########################################
# Start plotting
fig, ax = plt.subplots(figsize = (10,10))


########################
# Plot the Histogram   
hist, bin_edge = np.histogram(x_vals, weights = Weights, bins = Bins)
ax.bar(center_bins, hist/binwidth, width= binwidth, color='white', fill=False, edgecolor=colors,lw = 1.) 

########################
# Add the KDE
kernel = stats.gaussian_kde(x_vals, bw_method=0.1, weights=Weights)
# plot the KDE at these x-ticks
x_KDE = np.arange(0,40,0.5) 
center_KDEbins = (x_KDE[:-1] + x_KDE[1:])/2.
# Re-normalise the KDE
y_vals = kernel(center_KDEbins)*sum(hist)

# Plot it and fill the area underneath the curve :)
ax.plot(center_KDEbins, y_vals, lw=5, color=colors)
ax.fill_between(center_KDEbins, y1=0, y2=y_vals, color=colors, alpha = 0.1)

########################
# Bootstrap to obtain sampling errors
if bootstrap:
    indices = np.arange(len(x_vals))
    hist_vals = np.zeros((bootstraps, len(x_KDE)))  #center_bins
    for b in range(len(hist_vals)):
        print('boostrap b=', b)
        boot_index = np.random.choice(indices, size=len(indices), replace=True)

        kernel       = stats.gaussian_kde(x_vals[boot_index], bw_method=0.1, weights=Weights[boot_index])
        Hist, _      = np.histogram(x_vals[boot_index], bins=Bins,weights=Weights[boot_index],density=False)

        hist_vals[b] = kernel(x_KDE)*sum(Hist)

    # calculate 1- and 2- sigma percentiles
    percentiles = np.percentile(hist_vals, [15.89, 84.1, 2.27, 97.725], axis=0)
    median = np.percentile(hist_vals, [50], axis=0)
    
    ax.fill_between(x_KDE, percentiles[0],percentiles[1], alpha=0.4, color=colors, zorder = 11) # 1-sigma
    ax.fill_between(x_KDE, percentiles[2], percentiles[3],alpha=0.2, color=colors, zorder = 10) # 2-sgima    

    
#########################################
# plot values
ax.set_xlim(0,40)
ax.set_ylim(0,1)

ax.set_xlabel('$\mathcal{M}_c \ \mathrm{[M}_{\odot} \mathrm{]}$', fontsize = 30)
ax.set_ylabel('$\\frac{dN}{dt d\mathcal{M}_{c}}$ \ $\mathrm{[yr}^{-1}\mathrm{M}_{\odot}^{-1}\mathrm{]}$', fontsize = 30)

plt.savefig(fileLoc+'/MchirpO1detection_.pdf', bbox_inches='tight')

plt.show()<jupyter_output>boostrap b= 0
boostrap b= 1
boostrap b= 2
boostrap b= 3
boostrap b= 4
boostrap b= 5
boostrap b= 6
boostrap b= 7
boostrap b= 8
boostrap b= 9
boostrap b= 10
boostrap b= 11
boostrap b= 12
boostrap b= 13
boostrap b= 14
boostrap b= 15
boostrap b= 16
boostrap b= 17
boostrap b= 18
boostrap b= 19
boostrap b= 20
boostrap b= 21
boostrap b= 22
boostrap b= 23
boostrap b= 24
boostrap b= 25
boostrap b= 26
boostrap b= 27
boostrap b= 28
boostrap b= 29
boostrap b= 30
boostrap b= 31
boostrap b= 32
boostrap b= 33
boostrap b= 34
boostrap b= 35
boostrap b= 36
boostrap b= 37
boostrap b= 38
boostrap b= 39
boostrap b= 40
boostrap b= 41
boostrap b= 42
boostrap b= 43
boostrap b= 44
boostrap b= 45
boostrap b= 46
boostrap b= 47
boostrap b= 48
boostrap b= 49
{'fraction_non_alphanumeric': 0.17800498099357714, 'fraction_numerical': 0.04024118495215625, 'mean_word_length': 4.815548780487805, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '652084', 'license': 'MIT', 'n_tokens_mistral': 2943, 'n_tokens_neox': 2660, 'n_words': 721, 'path': 'examples/methods_paper_plots/chirpmass_distribution/chirpmass_dist.ipynb', 'repo_name': 'FloorBroekgaarden/COMPAS-2', 'size': '7629'}