Remnant-fits.md

@@ -4,6 +4,48 @@ To calculate the non-evolved remnant fits, we call the functions in `lalinferenc
 
 The conversions are handled in by the [`pesummary.gw.file.conversions` module](https://git.ligo.org/lscsoft/pesummary/-/blob/master/pesummary/gw/file/conversions.py) and specifically the [`pesummary.gw.file.conversion._Conversion` class](https://git.ligo.org/lscsoft/pesummary/-/blob/master/pesummary/gw/file/conversions.py#L534).
 
+To check the samples produced by PESummary are correct, we compared them to cbcBayesPostProc. We copied the following files from HAWK to CIT:
+
+```bash
+$ rsync -avz --rsh="gsissh" ligo.gravity.cf.ac.uk:/home/sebastian.khan/public_html/LVC/projects/o3/s190521g/run2-8s/1242442967.4472656-333631/H1L1V1/posterior_samples.dat .
+
+$ rsync -avz --rsh="gsissh" ligo.gravity.cf.ac.uk:/home/sebastian.khan/projects/o3/s190521g/paper-runs/run2-8s/run/posterior_samples/posterior_H1L1V1_1242442967.4472656-333631.hdf5 .
+```
+
+We then ran the following script:
+
+```python
+from pesummary.gw.file.read import read
+from pesummary.gw.file.standard_names import lalinference_map
+import numpy as np
+
+lalinf_samples = np.genfromtxt("/home/charlie.hoy/projects/pesummary_review/remnant_fits/non_evolved_average_precessing/posterior_samples.dat", names=True)
+PESUMMARY = read("/home/charlie.hoy/public_html/projects/pesummary_review/remnant_fits/non_evolved_average_precessing/comparison_to_master/samples/PhenomPv3HM_pesummary.dat")
+pesummary_samples = PESUMMARY.samples_dict
+
+reverse_map = {item: key for key, item in lalinference_map.items()}
+pesummary_dict = {}
+for param in pesummary_samples.keys():
+    if param in reverse_map.keys():
+        pesummary_dict[reverse_map[param]] = pesummary_samples[param]
+    else:
+        pesummary_dict[param] = pesummary_samples[param]
+
+# Check posterior samples
+not_included = []
+for param in ["l_peak_nonevol", "af_nonevol", "mf_nonevol"]:
+    if param in pesummary_dict.keys():
+        comparison = np.round(pesummary_dict[param], 8) == np.round(lalinf_samples[param], 8)
+        assert all(i == True for i in comparison)
+    else:
+        not_included.append(param)
+
+if len(not_included):
+    print("--------------------")
+    print("Parameters not included in this analysis: {}".format(", ".join(not_included)))
+    print("--------------------")
+```
+
 # Review comments
 
 ## David Keitel 20200304