Remnant-fits.md

@@ -181,7 +181,133 @@ Comments on the revised version
 
 ## example outputs
 
-🚧
+* https://ldas-jobs.ligo.caltech.edu/~charlie.hoy/projects/pesummary_review/remnant_fits/evolved_average_precessing/comparison_to_master/html/PhenomPv3HM_PhenomPv3HM_final_mass_evolved.html
+* https://ldas-jobs.ligo.caltech.edu/~charlie.hoy/projects/pesummary_review/remnant_fits/evolved_average_precessing/check_non_evolved_fits/comparison_to_master/home.html
+
+## Consistent non-evolved quantities
+
+* As the code base has changed to allow for the evaluation of the evolved spins, we check that we still find good agreement when we calculate non-evolved spins. 
+
+When we 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/evolved_average_precessing/check_non_evolved_fits/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", "e_rad_nonevol"]:
+    if param in pesummary_dict.keys() and param != "e_rad_nonevol":
+        comparison = np.round(pesummary_dict[param], 8) == np.round(lalinf_samples[param], 8)
+        assert all(i == True for i in comparison)
+        difference = pesummary_dict[param] - lalinf_samples[param]
+        print(
+            "Maximum difference for {}: {}".format(
+                param, np.max(np.abs(difference))
+            )
+        )
+    elif param == "e_rad_nonevol":
+        corrected_lalinf = lalinf_samples[param] * (pesummary_dict["m1_source"] / lalinf_samples["m1_source"])
+        comparison = np.round(pesummary_dict[param], 8) == np.round(corrected_lalinf, 8)
+        assert all(i == True for i in comparison)
+        difference = pesummary_dict[param] - corrected_lalinf
+        print(
+            "Maximum difference for {}: {}".format(
+                param, np.max(np.abs(difference))
+            )
+        )
+    else:
+        not_included.append(param)
+
+if len(not_included):
+    print("--------------------")
+    print("Parameters not included in this analysis: {}".format(", ".join(not_included)))
+    print("--------------------")
+```
+
+We get the following output:
+
+```bash
+$ python compare_output.py
+Maximum difference for l_peak_nonevol: 1.0045297926808416e-12
+Maximum difference for af_nonevol: 1.3800072196090696e-13
+Maximum difference for mf_nonevol: 7.958078640513122e-13
+Maximum difference for e_rad_nonevol: 4.1033842990145786e-13
+```
+
+Which is the same level of agreement as before. 
+
+## cbcBayesPostProc comparison
+
+* 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 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("IMRPhenomPv3HM_final_state_nrutils_evol.dat", names=True)
+PESUMMARY = read("/home/charlie.hoy/public_html/projects/pesummary_review/remnant_fits/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 = {}
+reverse_map["peak_luminosity_evolved"] = "l_peak_evol_avg"
+reverse_map["final_mass_evolved"] = "mf_evol_avg"
+reverse_map["tilt_1_evolved"] = "tilt1_isco"
+reverse_map["tilt_2_evolved"] = "tilt2_isco"
+reverse_map["phi_12_evolved"] = "phi12_isco"
+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 ["tilt1_isco", "tilt2_isco", "phi12_isco", "l_peak_evol_avg", "mf_evol_avg", "e_rad_evol_avg", "af_evol_avg"]:
+    if param in pesummary_dict.keys() and param != "e_rad_evol_avg":
+        comparison = np.round(pesummary_dict[param], 8) == np.round(lalinf_samples[param], 8)
+        difference = np.abs(pesummary_dict[param] - lalinf_samples[param])
+        """
+        print(
+            "Maximum difference for {}: {}".format(
+                param, np.max(np.abs(difference))
+            )
+        )
+```
+
+With output:
+
+```bash
+Maximum difference for tilt1_isco: 0.00016935480523350144
+Maximum difference for tilt2_isco: 0.00032114999367438557
+Maximum difference for phi12_isco: 0.0008060677213324041
+Maximum difference for l_peak_evol_avg: 1.0976890889136115e-05
+Maximum difference for mf_evol_avg: 2.745710679619151e-05
+Maximum difference for af_evol_avg: 9.770286555355057e-07
+```
 
 ## Review comments