SummaryReview-script.md

@@ -10,7 +10,7 @@ Can be seen here.
 
 # Review comments
 #### Step 1 (Run)
-We run `pesummary` and `cbcBayesPostProc` on the posterior file `.hdf5` indicated at [PESummary-Review#Notes](https://git.ligo.org/lscsoft/pesummary/-/wikis/PESummary-Review#notes)
+I run `pesummary` and `cbcBayesPostProc` on the posterior file `.hdf5` indicated at [PESummary-Review#Notes](https://git.ligo.org/lscsoft/pesummary/-/wikis/PESummary-Review#notes)
 
 ```bash
 summaryreview -w /home/nicola.delillo/public_html/pesummary-review/summaryreview -s /home/nicola.delillo/pesummary-review/samples/PROD0_posterior_samples.hdf5
@@ -27,7 +27,7 @@ $ lalsim.GetApproximantFromString('IMRPhenomPv2_NRtidal')
 79
 ```
 
-We edit the file running
+I edit the file running
 
 ```bash
 f = h5py.File('PROD0_posterior_samples.hdf5', 'a')
@@ -48,11 +48,7 @@ The single summary pages can be found at:
 In `pesummary` the credible interval for 1D histograms is calculated in a way that includes the 90% of the probability: [median - 5% percentile, median + 95%], as you can see from [plot.py#L106](https://git.ligo.org/lscsoft/pesummary/-/blob/master/pesummary/core/plots/plot.py#L106) and [plot.py#L107](https://git.ligo.org/lscsoft/pesummary/-/blob/master/pesummary/core/plots/plot.py#L107)
 in the source code.
 
-To review the credible interval produced, we use
-
-
-We have run this script that computes prints the credible intervals as output
-of pesummary attribute `.confidence_interval`
+To review the credible interval produced, have run this script that computes and  prints the credible intervals as output of pesummary attribute `.confidence_interval`
 
 
 ```python
@@ -116,13 +112,61 @@ L1_matched_filter_snr_angle [-0.07  0.06]
 a_1 [0.   0.11]
 ```
 
-These numbers have been compared to the reviewed results for GW190814 `combinedPHM` at the bottom  of this page: [https://git.ligo.org/publications/gw190814/gw190814-discovery/-/wikis/PE-result-table](https://git.ligo.org/publications/gw190814/gw190814-discovery/-/wikis/PE-result-table) 
+These numbers have been compared to the reviewed results for GW190814 `combinedPHM` at the bottom  of this page [1]
+
+
+* upper value `mass ratio`
+* lower value `total_mass_source`
+* lower value `mass_1_source`
+* lower limit `theta_jn `
+* both limits for `L1_matched_filter_snr_abs ` (0.05 difference for the lower limit)
+* Same for `H1_matched_filter_snr_abs `
+* Same for `V1_matched_filter_snr_abs` , with a 0.13 difference.
+
+
+My guess is that this is due to different versions of `np.percentile` (pesummary uses ). To cross-check this, I have personally run the script at the page linked [1]:
+
+```
+
+import numpy as np
+from pesummary.gw.file.read import read
+
+parameter_dict = {
+        "chirp_mass_source": "Chirp mass $\mathcal{M}/M_{\odot}$",
+        "total_mass_source": "Total mass $M/M_{\odot}$",
+        "mass_ratio": "Mass ratio $q$",
+        "mass_1_source": "Primary mass $m_{1}/M_{\odot}$",
+        "mass_2_source": "Secondary mass ${m_{2}/M_{\odot}}$",
+        "chi_eff": "Effective inspiral spin parameter $\chi_{eff}$",
+        "chi_p": "Effective precession parameter $\chi_{p}$",
+        "a_1": "Dimensionless primary spin magnitude $a_{1}$",
+        "luminosity_distance": "Luminosity Distance $D_{L}/\text{Mpc}$",
+        "redshift": "Source redshift $z$",
+        "theta_jn": "Inclination angle $\\theta_{\mathrm{JN}}$",
+        "L1_matched_filter_snr": "Signal to Noise ratio in \\llo $\\rho_{\mathrm{L}}$",
+        "H1_matched_filter_snr": "Signal to Noise ratio in \\lho $\\rho_{\mathrm{H}}$",
+        "V1_matched_filter_snr": "Signal to Noise ratio in \\virgo $\\rho_{\mathrm{V}}$",
+        "network_matched_filter_snr": "Network Signal to Noise ratio $\\rho_{\mathrm{HLV}}$",
+    }
+
+interested = ["PhenomPv3HM", "SEOBNRv4PHM", "combinedPHM"]
+for waveform_model in interested:
+    print("---------- {} ----------\n".format(waveform_model))
+    f = read(result_file_map[waveform_model])
+    f.generate_all_posterior_samples()
+    params = parameter_dict.keys()
+
+    for param in params:
+        samples = f.samples_dict[param]
+        print(
+            "{} = {}^{}_{}".format(
+                param, np.round(np.median(samples), 2),
+                np.round(np.percentile(samples, 95) - np.median(samples), 2),
+                np.round(np.median(samples) - np.percentile(samples, 5), 2)
+            )
+        )
+
+```
 
 
-* Check upper value `mass ratio`
-* Check lower value `total_mass_source`
+[1]: [https://git.ligo.org/publications/gw190814/gw190814-discovery/-/wikis/PE-result-table](https://git.ligo.org/publications/gw190814/gw190814-discovery/-/wikis/PE-result-table) . There are sum difference in the following parameters
\ No newline at end of file
-* Check lower value `mass_1_source`
-* Check lower limit `theta_jn `
-* Check both limits for `L1_matched_filter_snr_abs ` (0.05 difference for the lower limit)
-* Same issue with `H1_matched_filter_snr_abs `
-* Same issue with `V1_matched_filter_snr_abs` , with a 0.13 difference.