{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "expressed-speed",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib as mpl\n",
"import matplotlib.pyplot as plt\n",
"from pycbc.inference import io, models\n",
"from pycbc import waveform\n",
"from pycbc.detector import Detector\n",
"from pycbc import cosmology\n",
"\n",
"# PLOTTING OPTIONS\n",
"fig_width_pt = 3*246.0 # Get this from LaTeX using \\showthe\\columnwidth\n",
"inches_per_pt = 1.0/72.27 # Convert pt to inch\n",
"golden_mean = (np.sqrt(5)-1.0)/2.0 # Aesthetic ratio\n",
"fig_width = fig_width_pt*inches_per_pt # width in inches\n",
"fig_height = fig_width*golden_mean # height in inches\n",
"fig_size = [fig_width,fig_height]\n",
"\n",
"params = { 'axes.labelsize': 24,\n",
" 'font.family': 'serif',\n",
" 'font.serif': 'Computer Modern Raman',\n",
" 'font.size': 24,\n",
" 'legend.fontsize': 20,\n",
" 'xtick.labelsize': 24,\n",
" 'ytick.labelsize': 24,\n",
" 'axes.grid' : True,\n",
" 'text.usetex': True,\n",
" 'savefig.dpi' : 100,\n",
" 'lines.markersize' : 14,\n",
" 'figure.figsize': fig_size}\n",
"\n",
"mpl.rcParams.update(params)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "settled-motel",
"metadata": {},
"outputs": [],
"source": [
"def fn_to_waveform(fn): \n",
" # read in the raw posterior file and extract the max loglikelihood parameter value\n",
" with io.loadfile(fn, 'r') as fp: \n",
" data = fp.read_data() \n",
" psds = fp.read_psds() \n",
" cp = fp.read_config_file()\n",
" samples = fp.read_samples(list(fp['samples'].keys()))\n",
" idx = samples['loglikelihood'].argmax()\n",
" params = {p: samples[p][idx] for p in samples.fieldnames}\n",
" redshift = cosmology.redshift(params['distance'])\n",
" params['comoving_volume'] = cosmology.cosmological_quantity_from_redshift(redshift, 'comoving_volume')\n",
" # set up the model\n",
" model = models.read_from_config(cp, data=data, psds=psds)\n",
" model.update(**params)\n",
" _ = model.loglikelihood\n",
" pol = model.current_stats['maxl_polarization']\n",
" #compute the max loglikelihood waveform\n",
" wfs = model.waveform_generator.generate(polarization=pol,**model.current_params)\n",
" \n",
" whdata = {}\n",
" for det in model.data:\n",
" d = model.data[det]\n",
" d = d.copy()\n",
" d *= model.weight[det] / (4*model.psds[det].delta_f)**0.5\n",
" whdata[det] = d.to_timeseries()\n",
" \n",
" waveforms = {}\n",
" whitened_waveforms = {}\n",
" for detname in wfs:\n",
" hs = {}\n",
" whitend_hs = {}\n",
" det = Detector(detname)\n",
" fp, fc = det.antenna_pattern(model.current_params['ra'],\n",
" model.current_params['dec'],\n",
" pol,\n",
" model.current_params['tc'])\n",
" hp, hc = wfs[detname]\n",
" h = fp*hp + fc*hc\n",
" h.resize(len(model.weight[detname]))\n",
" wh = h.copy()\n",
" wh *= model.weight[detname] / (4*model.psds[detname].delta_f)**0.5\n",
" whitened_waveforms[detname] = wh.to_timeseries() \n",
" return whdata, whitened_waveforms,model"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "8691a3c5",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING:root:No sampling_params section read from config file\n",
"WARNING:root:WARNING: The following args are not being used for waveform generation: redshift, srcmchirp, comoving_volume, delta_tc, parity_mpvinverse, q, baseapprox, trigger_time\n",
"WARNING: fval is not bracketed by func(zmin)=3.630907107506232e-13 Mpc3 and func(zmax)=164042102234.23068 Mpc3. This means either there is no solution, or that there is more than one solution between zmin and zmax satisfying fval = func(z). [astropy.cosmology.funcs]\n",
"WARNING:astropy:fval is not bracketed by func(zmin)=3.630907107506232e-13 Mpc3 and func(zmax)=164042102234.23068 Mpc3. This means either there is no solution, or that there is more than one solution between zmin and zmax satisfying fval = func(z).\n",
"WARNING: fval is not bracketed by func(zmin)=3.630907107506232e-13 Mpc3 and func(zmax)=3763143105730.4414 Mpc3. This means either there is no solution, or that there is more than one solution between zmin and zmax satisfying fval = func(z). [astropy.cosmology.funcs]\n",
"WARNING:astropy:fval is not bracketed by func(zmin)=3.630907107506232e-13 Mpc3 and func(zmax)=3763143105730.4414 Mpc3. This means either there is no solution, or that there is more than one solution between zmin and zmax satisfying fval = func(z).\n",
"WARNING:root:No sampling_params section read from config file\n",
"WARNING:root:WARNING: The following args are not being used for waveform generation: redshift, srcmchirp, comoving_volume, delta_tc, q, trigger_time\n"
]
}
],
"source": [
"parity_path = '/work/yifan.wang/testingGR-o3b/snr/GW190521/convert.hdf'\n",
"parity_whdata, parity_whwaveform, parity_model = fn_to_waveform(parity_path)\n",
"\n",
"gr_path = '/work/yifan.wang/gw190521_xrerun/t1-mpi-gr/mpirun_output/config_files/convert.hdf'\n",
"_, gr_whwaveform,gr_model = fn_to_waveform(gr_path)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "d27e743c",
"metadata": {},
"outputs": [],
"source": [
"xlim=(-0.1, 0.1)\n",
"fig = plt.figure(figsize=(10.211706102117061, 6.311181454233049*1.6))\n",
"fig.subplots_adjust(hspace=.2)\n",
"ax = fig.add_subplot(311);\n",
"\n",
"tc = parity_model.static_params['trigger_time']\n",
"det = list(gr_whwaveform.keys())\n",
"wdata = parity_whdata[det[0]]\n",
"ax.plot(wdata.sample_times-tc, wdata, c='k', alpha=0.2, label='whitened data')\n",
"\n",
"wwf = parity_whwaveform[det[0]]\n",
"ax.plot(wwf.sample_times-tc, wwf,label='birefringence')\n",
"ax.plot(gr_whwaveform[det[0]].sample_times-tc, gr_whwaveform[det[0]],label='GR')\n",
"\n",
"ax.grid(True)\n",
"ax.set_xticklabels([])\n",
"ax.set_xlim(*xlim)\n",
"ax.set_ylim(-100,100)\n",
"ax.text(0.054, 70, 'LIGO Hanford', style='italic',fontsize=20,\n",
" bbox={'facecolor':'white', 'alpha':0.5, 'pad':5,'edgecolor':'none'})\n",
"ax.set_title('GW190521')\n",
"################\n",
"bx = fig.add_subplot(312);\n",
"wdata = parity_whdata[det[1]]\n",
"bx.plot(wdata.sample_times-tc, wdata, c='k', alpha=0.2, label='whitened data')\n",
" \n",
"wwf = parity_whwaveform[det[1]]\n",
"bx.plot(wwf.sample_times-tc, wwf,label='birefringence')\n",
"bx.plot(gr_whwaveform[det[1]].sample_times-tc, gr_whwaveform[det[1]],label='GR')\n",
"\n",
"bx.set_ylabel('Whitened Strain')\n",
"bx.set_xlim(*xlim)\n",
"bx.set_ylim(-100,100)\n",
"bx.text(0.0515, 70, 'LIGO Livingston', style='italic',fontsize=20,\n",
" bbox={'facecolor':'white', 'alpha':0.5, 'pad':0,'edgecolor':'none'})\n",
"bx.grid(True)\n",
"bx.set_xticklabels([])\n",
"\n",
"################\n",
"cx = fig.add_subplot(313);\n",
"wdata = parity_whdata[det[2]]\n",
"cx.plot(wdata.sample_times-tc, wdata, c='k', alpha=0.2, label='whitened data')\n",
" \n",
"wwf = parity_whwaveform[det[2]]\n",
"cx.plot(wwf.sample_times-tc, wwf,label='birefringence')\n",
"cx.plot(gr_whwaveform[det[2]].sample_times-tc, gr_whwaveform[det[2]],label='GR')\n",
"cx.text(0.07, 70, 'Virgo', style='italic',fontsize=20,\n",
" bbox={'facecolor':'white', 'alpha':0.5, 'pad':5,'edgecolor':'none'})\n",
"\n",
"cx.set_xlim(*xlim)\n",
"cx.set_ylim(-100,100)\n",
"cx.set_xlabel('Time - '+str(tc)+' [s]')\n",
"cx.legend(loc='lower right', ncol=3,framealpha=0)\n",
"fig.savefig('0521_gw190521.pdf',bbox_inches='tight')"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "a2769cdb",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'inclination': 1.3773388232444717,\n",
" 'maxl_polarization': 5.132211422080623,\n",
" 'srcmchirp': 59.988695836431965,\n",
" 'logwt': -28046.517628006437,\n",
" 'V1_optimal_snrsq': 4.645117051883857,\n",
" 'parity_mpvinverse': 571.5796601060491,\n",
" 'distance': 1333.6349776943568,\n",
" 'spin1_a': 0.9116850950619036,\n",
" 'dec': 0.5513436858058255,\n",
" 'delta_tc': 0.0061986454908000005,\n",
" 'spin1_polar': 2.076365152215487,\n",
" 'spin2_azimuthal': 6.273223231326342,\n",
" 'spin2_polar': 0.9148528363767893,\n",
" 'ra': 3.3490655420980504,\n",
" 'logprior': -49.260229458075116,\n",
" 'spin2_a': 0.09566031928959656,\n",
" 'coa_phase': 2.930192323180816,\n",
" 'loglr': 122.35747765652349,\n",
" 'L1_optimal_snrsq': 163.8200365758998,\n",
" 'loglikelihood': -28001.80399824368,\n",
" 'logjacobian': 0.0,\n",
" 'q': 4.542030730356726,\n",
" 'maxl_loglr': 125.90973288111286,\n",
" 'spin1_azimuthal': 2.6322431174435468,\n",
" 'H1_optimal_snrsq': 69.48418153589377,\n",
" 'comoving_volume': 5018380778.811694,\n",
" 'approximant': 'mpvnosmall',\n",
" 'baseapprox': 'IMRPhenomXPHM',\n",
" 'f_lower': 20.0,\n",
" 'f_ref': 20.0,\n",
" 'trigger_time': 1242442967.4473,\n",
" 'spin1x': -0.6963825154031204,\n",
" 'spin1y': 0.3889331472174398,\n",
" 'spin1z': -0.44153381727155666,\n",
" 'spin2x': 0.07580444315196702,\n",
" 'spin2y': -0.0007551945954393175,\n",
" 'spin2z': 0.05834391798748846,\n",
" 'tc': 1242442967.4534986,\n",
" 'redshift': 0.2556724782843754,\n",
" 'mass1': 194.34779119906074,\n",
" 'mass2': 42.78874422846471}"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"parity_model.current_params"
]
},
{
"cell_type": "markdown",
"id": "7476300f",
"metadata": {},
"source": [
"# SNR"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "6bd3b807",
"metadata": {},
"outputs": [],
"source": [
"import h5py\n",
"from pycbc.conversions import snr_from_loglr\n",
"nongr = h5py.File(parity_path,'r')\n",
"gr = h5py.File(gr_path,'r')"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "2b6eb9a3",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"15.643367774013592"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"snr_from_loglr(np.max(nongr['samples']['loglr']))"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "358ab6ca",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"15.312425524512493"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"snr_from_loglr(np.max(gr['samples']['loglr']))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "eae483fd",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}