Merge pull request #77 from legend-exp/refactor

small docu additions

docs/src/api.md

@@ -1,10 +1,19 @@
 # API
 
+Table of contents:
+
+```@contents
+Pages = ["api.md"]
+Depth = 3
+```
+
+## Functions
 ```@index
 Pages = ["api.md"]
 Order = [:function]
 ```
 
+## Documentation
 ### analysis.jl
 ```@docs
 ZeroNuFit.Analysis.retrieve_real_fit_results

docs/src/likelihood.md

@@ -22,7 +22,7 @@ In case no events are found in a given partition _k_, the above Likelihood expre
 
 ```math
 \begin{aligned}
-    \mathcal{L}(\Gamma) = \prod_k \textrm{Pois}(s_k+b_k) 
+    \mathcal{L}(\Gamma,\, \boldsymbol{BI},\,\boldsymbol{\theta}|D) = \prod_k \textrm{Pois}(s_k+b_k) 
 \end{aligned}
 ```
 
@@ -142,6 +142,58 @@ The above products, then, can be expressed again as
 \end{aligned}
 ```
 
+### Alternative background shapes
+
+The default background modeling shape is a flat function:
+
+```math
+\begin{aligned}
+  f_{\rm flat}(E) = b_{\rm k} \cdot p_{\rm b,\,flat}(E) = BI \cdot \Delta E \cdot \mathcal{E}_{\rm k} \cdot \underbrace{\frac{1}{K_{\rm flat}}}_{p_{\rm b,\,flat}(E)}
+\end{aligned}
+```
+
+where $E$ is the energy and $K_{\rm flat}$ is a normalization factor that accounts for the net width of the fit window, accounting for any removed gap within it.
+
+In the linear case, we can model the background as
+
+```math
+\begin{aligned}
+  f_{\rm lin}(E) = b_{\rm k} \cdot p_{\rm b,\,lin}(E) = BI \cdot \Delta E \cdot \mathcal{E}_{\rm k} \cdot \underbrace{\left(1+ \frac{m_{\rm lin} \cdot (E-E_{\rm 0})}{E_{\rm max}-E_{\rm 0}} \right)\cdot \frac{1}{K_{\rm lin}}}_{p_{\rm b,\,lin}(E)}
+\end{aligned}
+```
+
+where $m_{\rm lin}$ is the slope of the linear function, $E_{\rm 0}$ ($E_{\rm max}$) is the starting (ending) energy value of the fit window, and $K_{\rm lin}$ is the normalization factor.
+
+In the exponential case, we can model the background as
+
+```math
+\begin{aligned}
+  f_{\rm exp}(E) = b_{\rm k} \cdot p_{\rm b,\,exp}(E) = BI \cdot \Delta E \cdot \mathcal{E}_{\rm k} \cdot \underbrace{e ^ {\left(E-E_{\rm 0}\right)\cdot \frac{m_{\rm exp}}{\Delta E} } \cdot \frac{1}{K_{\rm exp}}}_{p_{\rm b,\,exp}(E)} 
+\end{aligned}
+```
+
+with corresponding slope $m_{\rm exp}$ and normalization factor $K_{\rm exp}$.
+
+The normalization factors can be expressed in a general form in the following way:
+
+```math
+\begin{aligned}
+    \begin{cases}
+      K_{\rm flat} = \Delta E\\[15pt]
+      K_{\rm lin} = \Delta E \cdot \left(1 - \frac{m_{\rm lin}\cdot E_{\rm 0}}{E_{\rm max}-E_{\rm 0}} \right) + m_{\rm lin}\cdot \frac{\sum_{i} \left( E_{\rm h,\,i}^2 - E_{\rm l,\,i}^2 \right)}{2\left(E_{\rm max}-E_{\rm 0} \right )}\\[15pt]
+      K_{\rm exp} = \frac{E_{\rm max}-E_{\rm 0}}{m_{\rm exp}} \cdot \left[\sum_{\rm i}e^{(E_{\rm h,i}-E_{\rm 0})\cdot \frac{m_{\rm exp}}{E_{\rm max}-E_{\rm 0}}} - \sum_{\rm i}e^{(E_{\rm l,i}-E_{\rm 0})\cdot \frac{m_{\rm exp}}{E_{\rm max}-E_{\rm 0}}} \right]
+    \end{cases}  
+\end{aligned}
+```
+
+where $E_{\rm l,i}$ ($E_{\rm h,i}$) is the starting (ending) energy value of the sub-windows defined within the analysis fit window, accounting for any excluded gap in the fit window. 
+If no gaps are present, then $E_{\rm l,i}\equiv E_{\rm 0}$ and $E_{\rm h,i}\equiv E_{\rm max}$.
+In particular, $\Delta E= \sum_{\rm i} \left(E_{\rm h,i}-E_{\rm l,i} \right)$.
+
+Notice that $K_{\rm lin} \rightarrow K_{\rm flat}$ and $K_{\rm exp} \rightarrow K_{\rm flat}$ in the limit of $m_{\rm lin} \rightarrow 0$ and $m_{\rm exp} \rightarrow 0$, respectively. 
+The free parameters in the linear and exponential case are $m_{\rm lin}$ and $m_{\rm exp}$ (other than the background indices, $BI$).
+For all these parameters, we used a uniform prior as default, while prior ranges can be modified by the user via the configuration file. 
+
 ### Posterior distributions and marginalization
 
 The combined posterior probability density function is calculated according to Bayes’ theorem as:

docs/src/toys.md

@@ -7,46 +7,60 @@ Pages = ["toys.md"]
 Depth = 3
 ```
 
-## Running toys for posterior predictive distribution studies
+## Posterior predictive distribution studies
 
 In order to check for any mis-modeling in the fits and to further understand the derived half-life limits, we performed "sensitivity" studies. 
 In a Bayesian context, the sensitivity is related to the concept of posterior predictive distributions.
 The prior predictive distribution is the expected distribution of data coming from a future experiment identical to that performed and repeated under the same conditions.
 This marginalizes the uncertainty in the model parameters, based on their prior distributions. 
 Alternatively, the posterior predictive distribution weights the data by the posterior obtained from the analysis.
-If the original data were modeled appropriately, then fake data generated under the given model should distribute similarly to the original data.
+If the original data were modeled appropriately, then fake data generated under the given model should distribute similarly to the original data (see _A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin, "Bayesian Data Analysis", 3rd ed. 798, Chapman & Hall, 2013_ for further details).
 
 Considering the observed signal counts as an observable of the data, we can thus extract "sensitivity curves" that can be interpreted as the distribution of expected future limits derived from repeating the identical experiment.
 This is distinct from a frequentist sensitivity since uncertainty on nuisance parameters are marginalized over.
 
+Sensitivity studies can be run over $N_{\rm tot}$ toy spectra generated by sampling the posterior marginalized distributions of nuisance parameters obtained in the no-signal hypothesis fit, i.e. in a fit where the signal contribution was excluded from the modeling of events. 
+In particular, starting from the posterior distributions derived for all nuisance parameters via the fit over original data, $p(\theta | D)$, data were sampled from a model's posterior distribution as follows:
+
+```math
+\begin{aligned}
+    p(y | D) = \int p(y | \theta) p(\theta | D) d\theta
+\end{aligned}
+```
+
+where $y$ are the generated data, $D$ are the observed data, and $\theta$ are the model parameters (but the signal). 
+This is achieved by drawing samples of $\theta$ from the posterior distribution $p(\theta | D)$, and subsequently generating datasets $y$ based on the likelihood model $p(y | \theta)$. 
+In particular, for each partition, the expected number of background events was calculated using Poisson sampling, $n_{\rm b} \sim \text{Pois}(\mu_b)$, so that $n_{\rm b}$ events could be uniformly generated in the considered fit energy window to build a toy energy spectrum.
+Each generated toy energy spectrum can be later fitted with a signal+background model.
+
+## Running toys
+
 A module `sensitivity.jl` is present for generating toys and running sensitivity studies. The script can be run as
 
 ```
 $ julia sensitivity.jl -c config_fake_data.json -i N
 ```
 
-where `N` is an integer number corresponding to the toy index.
-The command can be run in an external bash script for looping over this index.
-
+where $N=\{1,\,...,\,N_{\rm tot}\}$ is an integer number corresponding to the toy index.
 The input config file (`config_fake_data.json`) has the following entries:
 
 ```
 {
     "path_to_fit": "output/fit_gerda_phIandphII_NoSignal/",
     "best_fit": false,
+    "low_stat": true,
     "seed": null
 }
 
 ```
 
 where
 - `"path_to_fit"` is the path to the already performed fit over real data;
-- `"best_fit": true` if we want to fix the paramaters to the best fit;
+- `"best_fit": true` if we want to fix the paramaters to the best fit coming from the fit under the no-signal hypothesis;
+- `"low_stat": true` to run 5 MCMC chains with $10^5$ steps each; if `false`, then the statistics used for the original fit under the no-signal hypothesis is used;
 - `"seed": null` if we want a random seed when generating fake data, otherwise you can fix it to an Int value.
 
-Any information about the signal being included or not in the fit of real data, was saved and retrieved from the output JSON file with results.
-
-Below, we show an example of bash file used for running sensitivity studies as multiple jobs on NERSC:
+Below, we show an example of a bash script used for running sensitivity studies as multiple jobs on NERSC:
 
 ```bash
 #!/bin/bash                                                                                                                                                 
@@ -69,15 +83,15 @@ wait
 ```
 
 
-## Using already existing toys to test alternative models
-Another way to run the code is present if, for instance, an user wants to use toy data generated according to one model but fit them with another model.
-In this case, the path to the folder containing the already existing JSON files with toy data has to be provided together with the toy index:
+## Testing alternative models using already existing toys 
+Another way to run the code is present if, for instance, an user wants to use toy data generated according to one model and fit them with a different model.
+In this case, the path to the folder containing the already existing JSON files with toy energy events previously generated has to be provided together with the toy index:
 
 ```
 julia sensitivity.jl -c config_fake_data.json -i N --path_to_toys path_to_your_toys
 ```
 
-Below, an updated version of a bash file that can be used for retrieving multiple existing toy data and running sensitivity studies as multiple jobs on NERSC:
+Below, an example of a bash script that can be used for retrieving multiple existing toy data and running sensitivity studies as multiple jobs on NERSC:
 
 ```bash
 #!/bin/bash