add formula in index
@@ -52,6 +52,59 @@ Note that in Gammapy, 2D image analyses are done with 3D cubes with a single
energy bin, e.g. for modeling and fitting,
see the `2D map analysis tutorial <./tutorials/image_analysis.html>`__.
+To analyse multiple runs, you can either stack the datasets together, or perform
+a joint fit across multiple datasets.
+
+Stacking Multiple Datasets
+==========================
+
+Stacking datasets implies that the counts, background and reduced IRFs from all the
+runs are binned together to get one final dataset for which a likelihood is
+computed during the fit.
+
+Counts, background and exposure (lying outside the masked regions are simply summed),
+while the energy dispersion and point spread function are weighed by the exposure
+and then summed.
+
+For the model evaluation, an important factor that needs to be accounted for is
+that the energy threshold changes between obseravtions.
+To ensure that the npred (ie, the predicted number of counts) on the stacked
+dataset is the sum expected by stacking the npred of the individual runs,
+the exposure weighted IRFs are computed.
+The mask_safe from each dataset is applied on the respective reconstructed energy axis
+of the energy dispersion matrix, and the masked matrices are combined.
+Values lying outside the safe mask of each dataset are lost.
+
+The stacking of 2 datasets is implemented as follows.
+Here, :math:`k` denotes a bin in reconstructed energy and
+:math:`j = {1,2}` is the dataset number
+
+
+================= ================================== ==================================================================================
+Dataset attribute Behaviour Implementation
+================= ================================== ===================================================================================
+`livetime` Sum of individual livetimes \overline{t} = \sum_j t_j
+`mask_safe` Pixels added with `OR` operation \overline{\epsilon_k} = \wedge_j \epsilon_{jk}
+`mask_fit` Dropped
+`counts` Summed outside exclusion regions \overline{\mathrm{n_{on}}_k} = \sum_j \mathrm{n_{on}}_{jk} \cdot \epsilon_{jk}
+`background` Summed outside exclusion regions \overline{\mathrm{bkg}_k} = \sum_j \mathrm{bkg}_{jk} \cdot \epsilon_{jk}
+`exposure` Summed outside exclusion regions \overline{\mathrm{aeff}}_l = \frac{\sum_{j}\mathrm{aeff}_{jl} \cdot t_j}{\overline{t}}
+`psf` Weighed by exposure and summed
+`edisp` Weighed by exposure, mask applied
+ on reco energy, and summed
+
+
+
+
+It is important to keep in mind that:
+- Stacking happens in-place, ie, dataset1.stack(dataset2) will overwrite dataset
+- To properly handle masks, it is necessary to stack onto an empty dataset.
+- Stacking only works for maps with equivalent geometry.
+Two geometries are called equivalent if one is exactly the same as,
+or can be obtained from a cutout of, the other.
+- A stacked analysis is reasonable only when adding runs taken by the same instrument.
+
+
Getting Started
===============
