draft: debug the redshift probability distribution

• Fix an error where the probability for the redshift distribution was always being reported as flat-in-comoving volume even when it was not
• Add normalization to the redshift part of pdraw
• Minor bugfix with PEsummary format

generate-mock-PE.py

@@ -24,16 +24,16 @@ import os
 # TODO: turn these settings into a config file
 
 # output label/path
-run_label = 'with_z_evo_PE_format_lalprior_69_evs_'
-pesummary_format=False
+run_label = 'with_z_evo_PE_format_lalprior_69_evs_z_debug'
+pesummary_format=True
 
 # how many events do we want in the end
-num_events = 1
-num_realizations = 1
+num_events = 69
+num_realizations = 20
 num_samples_per_event = 5000
 
 # do we want a corresponding injection set?
-create_injs=True
+create_injs=False
 
 # what parameterization would we like our samples in?
 # options are: "detectorframem1m2z", "detectorframeMcetarho",
@@ -167,7 +167,10 @@ def draw_true_from_plaw(mmin,mmax,alpha,beta,lambda_z,num_draws,
     zt = inv_zdist_O3(np.random.rand(num_draws))
 
     if return_draw_prob:
-        z_prob = cosmo_utils.ppop_z_func(cosmo_dict)(zt)
+        z_prob = cosmo_utils.ppop_z_func(cosmo_dict,lambda_z)(zt)
+        #zarr = np.linspace(0,zmax,num=1000)
+        #znorm = np.trapz(cosmo_utils.ppop_z_func(cosmo_dict,lambda_z)(zarr),zarr)
+        #z_prob /= znorm 
         draw_prob = powerlaw(m1t,-alpha,high=mmax,low=mmin) * powerlaw(q,beta,high=1.,low=mmin/m1t) * z_prob
         return m1t,m2t,zt, draw_prob
     else:
@@ -263,7 +266,7 @@ for j in tqdm(range(num_realizations)):
     if pesummary_format:
         if not os.path.exists(f'out/{run_label}'):
             os.makedirs(f'out/{run_label}')
-            os.makedirs(f'out/{run_label}/realization_{f}')
+            os.makedirs(f'out/{run_label}/realization_{j}')
         for event in ds:
             evnum = int(event[6:])
             sname = f"S2301{evnum:02}ab"
@@ -273,13 +276,14 @@ for j in tqdm(range(num_realizations)):
                 for param in ds.parameter.values:
                     group_pe.create_dataset(transforms.pesummary_key_map[param],
                                          data=ds[event].sel(parameter=param).values)
-                for param in true_vals_ds.parameter.values:
-                    group_true.create_dataset(transforms.pesummary_key_map[param],
-                                         data=true_vals_ds[event].sel(parameter=param).values)
+                if j<10:
+                    for param in true_vals_ds.parameter.values:
+                        group_true.create_dataset(transforms.pesummary_key_map[param],
+                                             data=true_vals_ds[event].sel(parameter=param).values)
 
 
 if create_injs:
-    m1i,m2i,zi,pdraw = draw_true_from_plaw(1,200,2.35, beta_q,
+    m1i,m2i,zi,pdraw = draw_true_from_plaw(mmin,200,2.35, beta_q,
                                            redshift_plaw,
                                            num_draws=int(2e7),
                                            return_draw_prob=True,
@@ -294,7 +298,7 @@ if create_injs:
 
     # put them in same format as Reed does
     # adaped from https://git.ligo.org/reed.essick/o1-o2-semianalytic-and-o3-real-injections/-/blob/main/create-multi-run-mixture#L159
-    with h5py.File(f"out/{run_label}_O3_L1_actual_injections.h5","w") as ff:
+    with h5py.File(f"out/{run_label}_mmin_mmin_no_norm_O3_L1_actual_injections.h5","w") as ff:
         group = ff.create_group("injections")
 
         group.create_dataset("mass1_source", data=m1i)