Commit e3a953307f1f9f1ea2a08bd96513a045fdef3fd4
1 parent
1ebf0992
Added tech note about RoI handling in stacked analysis (#1819)
Showing
1 changed file
with
205 additions
and
0 deletions
doc/dev/tn/gammalib_tn0002_roi_stacked.tex
0 → 100644
| 1 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 2 | +% Handling of Regions of Interest in stacked analysis | |
| 3 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 4 | + | |
| 5 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 6 | +% Definitions for manual package | |
| 7 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 8 | +\newcommand{\task}{\mbox{GammaLib}} | |
| 9 | +\newcommand{\this}{\mbox{\tt \task}} | |
| 10 | +\newcommand{\shorttype}{\mbox{}} | |
| 11 | +\newcommand{\doctype}{\mbox{RoI handling in stacked analysis}} | |
| 12 | +\newcommand{\reference}{\mbox{GammaLib-TN0002}} | |
| 13 | +\newcommand{\version}{\mbox{1.0}} | |
| 14 | +\newcommand{\calendar}{\mbox{20 July 2016}} | |
| 15 | +\newcommand{\auth}{\mbox{J\"urgen Kn\"odlseder}} | |
| 16 | +\newcommand{\approv}{\mbox{J\"urgen Kn\"odlseder}} | |
| 17 | + | |
| 18 | + | |
| 19 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 20 | +% Document definition | |
| 21 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 22 | +\documentclass{article}[12pt,a4] | |
| 23 | +\usepackage{epsfig} | |
| 24 | +\usepackage{technote} | |
| 25 | + | |
| 26 | + | |
| 27 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 28 | +% Begin of document body | |
| 29 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 30 | +\begin{document} | |
| 31 | +\frontpage | |
| 32 | + | |
| 33 | + | |
| 34 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 35 | +% Scope | |
| 36 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 37 | +\section{Scope} | |
| 38 | + | |
| 39 | +This technical note provides considerations for handling of Regions of Interest in a stacked | |
| 40 | +analysis. | |
| 41 | + | |
| 42 | + | |
| 43 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 44 | +% Formulating the problem | |
| 45 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 46 | +\section{Formulating the problem} | |
| 47 | + | |
| 48 | +Each analysis of an observation requires the definition of a Region of Interest (RoI) which delimits | |
| 49 | +the data space that is used for the analysis. | |
| 50 | +The RoI can be a spatial selection (for example a circle on the sky) or an energy selection, or a | |
| 51 | +combination of both. | |
| 52 | +Let's denote in the following the data space coordinates by $\vec{d}$. | |
| 53 | +The expected number of events in an observation $i$ after RoI selection is then given by | |
| 54 | +\begin{equation} | |
| 55 | +e_i(\vec{d}) = w_i(\vec{d}) \times M_i(\vec{d}) | |
| 56 | +\end{equation} | |
| 57 | +where $M_i(\vec{d})$ is the expected number of events (computed from a model) before the | |
| 58 | +RoI selection, and $w_i(\vec{d})$ is the RoI window function that is $1$ inside the RoI and $0$ | |
| 59 | +outside. | |
| 60 | +For stacked analysis the data space is represented by a counts cube, and the RoI can only partially | |
| 61 | +overlap with some of the cube bins. | |
| 62 | +The partial overlap can for example be spatially or in energy. | |
| 63 | +The window function can then be generalised into an function given the fractional overlap of a | |
| 64 | +counts cube bin with the RoI, and this function takes values between $0$ and $1$. | |
| 65 | + | |
| 66 | +Let's now consider the stacking of two observations into a single counts cube (the reasoning can | |
| 67 | +be easily extended to an arbitrary number of observations). | |
| 68 | +The expected number of events in the stacked counts cube is then given by | |
| 69 | +\begin{equation} | |
| 70 | +e(\vec{d}) = w_1(\vec{d}) M_1(\vec{d}) + w_2(\vec{d}) M_2(\vec{d}) \, . | |
| 71 | +\label{eq:added} | |
| 72 | +\end{equation} | |
| 73 | +We are now seeking for a weighting function $w(\vec{d})$ of the stacked observation so that | |
| 74 | +we can simply write the expected number of events as the product between the weighting | |
| 75 | +function and a stacked model, i.e. | |
| 76 | +\begin{equation} | |
| 77 | +e(\vec{d}) = w(\vec{d}) \times \left( M_1(\vec{d}) + M_2(\vec{d})\right) \, . | |
| 78 | +\label{eq:stacked} | |
| 79 | +\end{equation} | |
| 80 | +Combining Eq.~(\ref{eq:added}) and Eq.~(\ref{eq:stacked}) results in | |
| 81 | +\begin{eqnarray} | |
| 82 | +w(\vec{d}) \times \left( M_1(\vec{d}) + M_2(\vec{d})\right) & = & w_1(\vec{d}) M_1(\vec{d}) + w_2(\vec{d}) M_2(\vec{d}) \nonumber \\ | |
| 83 | +w(\vec{d}) & = & \frac{w_1(\vec{d}) M_1(\vec{d}) + w_2(\vec{d}) M_2(\vec{d})} {M_1(\vec{d}) + M_2(\vec{d})} \nonumber \\ | |
| 84 | +w(\vec{d}) & = & w_1(\vec{d}) \frac{M_1(\vec{d})}{M_1(\vec{d}) + M_2(\vec{d})} + | |
| 85 | +w_2(\vec{d}) \frac{M_2(\vec{d})}{M_1(\vec{d}) + M_2(\vec{d})} | |
| 86 | +\label{eq:weight} | |
| 87 | +\end{eqnarray} | |
| 88 | +which means that the weighting function for the stacked observation is the sum of the weighting | |
| 89 | +function of the individual observations, weighted by the relative contribution of the model in the | |
| 90 | +individual observations. | |
| 91 | +The obvious problem is that the relative contribution of the model in the individual observations | |
| 92 | +is in not available in a stacked analysis. | |
| 93 | + | |
| 94 | + | |
| 95 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 96 | +% Ignoring partial overlap | |
| 97 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 98 | +\section{Ignoring partial overlap} | |
| 99 | + | |
| 100 | +Suppose we can ignore the effects of partial overlap of the RoI with counts cube bins (we could | |
| 101 | +accomplish this in principle by redefining the bin boundaries so that there is not partial overlap | |
| 102 | +problem). | |
| 103 | +In that case the window functions of individual observations $w_i(\vec{d})$ take the values $0$ | |
| 104 | +and $1$. | |
| 105 | +We can then distinguish four cases for our example of stacking two observations: | |
| 106 | +\begin{equation} | |
| 107 | +w(\vec{d}) = \left \{ | |
| 108 | + \begin{array}{c l} | |
| 109 | + \displaystyle | |
| 110 | + 1 & w_1(\vec{d})=1, w_2(\vec{d})=1 \\ | |
| 111 | + \\ | |
| 112 | + \displaystyle | |
| 113 | + \frac{M_1(\vec{d})}{M_1(\vec{d})+M_2(\vec{d})} & w_1(\vec{d})=1, w_2(\vec{d})=0 \\ | |
| 114 | + \\ | |
| 115 | + \displaystyle | |
| 116 | + \frac{M_2(\vec{d})}{M_1(\vec{d})+M_2(\vec{d})} & w_1(\vec{d})=0, w_2(\vec{d})=1 \\ | |
| 117 | + \\ | |
| 118 | + \displaystyle | |
| 119 | + 0 & w_1(\vec{d})=0, w_2(\vec{d})=0 | |
| 120 | + \end{array} | |
| 121 | + \right . | |
| 122 | +\label{eq:cases} | |
| 123 | +\end{equation} | |
| 124 | +For all data space coordinates that fall in both RoIs the stacked weighting function is 1, while for all | |
| 125 | +coordinates that fall outside both RoIs the stacked weighting function is 0. | |
| 126 | +In case that the data space coordinate $\vec{d}$ falls in the RoI of the first observation but outside the | |
| 127 | +RoI of the second observation the weighting function becomes | |
| 128 | +\begin{equation} | |
| 129 | +w(\vec{d}) = \frac{M_1(\vec{d})}{M_1(\vec{d})+M_2(\vec{d})} | |
| 130 | +\end{equation} | |
| 131 | +which is smaller than 1 for all coordinates $\vec{d}$ for which $M_2(\vec{d}) > 0$. | |
| 132 | +These coordinates correspond to zones that fall {\bf outside the RoI of the second observation}, | |
| 133 | +hence the reduction of the weighting function is due to the spill over of events from the second | |
| 134 | +observation coming from outside its RoI into the stacked observation. | |
| 135 | + | |
| 136 | +A solution to this problem is to avoid the spill over, and to set $M_2(\vec{d}) = 0$ for all | |
| 137 | +$\vec{d}$ outside the RoI. | |
| 138 | +In other terms | |
| 139 | +\begin{equation} | |
| 140 | +\tilde{M_i}(\vec{d}) = w_i(\vec{d}) \times M_i(\vec{d}) | |
| 141 | +\end{equation} | |
| 142 | +The issue here is that this cannot be done in an exact way for sky models, since they are defined | |
| 143 | +in the physics space and not the data space, hence PSF blurring or energy redistribution cannot | |
| 144 | +be taken into account properly. | |
| 145 | +It nevertheless is probably still the best (maybe the only?) solution to limit the effect of the spill | |
| 146 | +over problem. | |
| 147 | + | |
| 148 | + | |
| 149 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 150 | +% Taking into account partial overlap | |
| 151 | +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 152 | +\section{Taking into account partial overlap} | |
| 153 | + | |
| 154 | +In case we want to take into account the partial overlap of counts cube bins with the RoIs we | |
| 155 | +need to compute Eq.~(\ref{eq:weight}). | |
| 156 | +In general, the model can be factorised into | |
| 157 | +\begin{equation} | |
| 158 | +M_i(\vec{d}) = R_i(\vec{d}) \times \tau_i | |
| 159 | +\end{equation} | |
| 160 | +where $R_i(\vec{d})$ is an event rate and $\tau_i$ is the livetime of the observation $i$. | |
| 161 | +Equation \ref{eq:weight} can then be rewritten as | |
| 162 | +\begin{equation} | |
| 163 | +w(\vec{d}) = | |
| 164 | +w_1(\vec{d}) \frac{R_1(\vec{d}) \times \tau_1} | |
| 165 | +{R_1(\vec{d}) \times \tau_1 + R_2(\vec{d}) \times \tau_2} + | |
| 166 | +w_2(\vec{d}) \frac{R_2(\vec{d}) \times \tau_2} | |
| 167 | +{R_1(\vec{d}) \times \tau_1 + R_2(\vec{d}) \times \tau_2} | |
| 168 | +\label{eq:overlap} | |
| 169 | +\end{equation} | |
| 170 | +We can still apply the logic of setting $R_i(\vec{d}) = 0$ outside the RoI, but there will be bins | |
| 171 | +(specifically in energy) for which the RoIs of both observations will have partial overlap, and we | |
| 172 | +need to compute Eq.~(\ref{eq:overlap}) explicitly. | |
| 173 | + | |
| 174 | +For this we may consider some special cases. | |
| 175 | +The first is that $R_1(\vec{d}) \approx R_2(\vec{d})$ which would be the case if both observations | |
| 176 | +share roughly the same pointing and have the same energy threshold. | |
| 177 | +We can then drop the $R_i(\vec{d})$ from Eq.~(\ref{eq:overlap}) and get | |
| 178 | +\begin{equation} | |
| 179 | +w(\vec{d}) \approx w_1(\vec{d}) \frac{\tau_1}{\tau_1 + \tau_2} + w_2(\vec{d}) \frac{\tau_2}{\tau_1 + \tau_2} | |
| 180 | +\end{equation} | |
| 181 | +hence the weighting function for the stacked analysis is the sum of the weighting functions of the | |
| 182 | +observations, weighted by their respective livetimes. | |
| 183 | + | |
| 184 | +Another special case is that there is negligible overlap between both observations, hence | |
| 185 | +$R_1(\vec{d}) \approx 0$ where $R_2(\vec{d})$ is non-zero, and vice versa. | |
| 186 | +Then all fractions in Eq.~(\ref{eq:overlap}) become unity and we get | |
| 187 | +\begin{equation} | |
| 188 | +w(\vec{d}) \approx w_1(\vec{d}) + w_2(\vec{d}) \, . | |
| 189 | +\end{equation} | |
| 190 | + | |
| 191 | +We can also consider the special case $R_1(\vec{d}) \times \tau_1 \gg R_2(\vec{d}) \times \tau_2$ | |
| 192 | +for which Eq.~(\ref{eq:overlap}) | |
| 193 | +is approximated by | |
| 194 | +\begin{equation} | |
| 195 | +w(\vec{d}) \approx w_1(\vec{d}) + w_2(\vec{d}) \frac{R_2(\vec{d}) \times \tau_2}{R_1(\vec{d}) \times \tau_1} | |
| 196 | +\approx w_1(\vec{d}) \, . | |
| 197 | +\end{equation} | |
| 198 | + | |
| 199 | +This illustrates that the weighting function taking into account partial overlaps can take different forms | |
| 200 | +dependent on the relative weights of the $R_i(\vec{d})$ or $M_i(\vec{d})$, and there is no obvious way | |
| 201 | +of how to do this properly without knowing the individual $M_i(\vec{d})$ or making an assumption about | |
| 202 | +them. | |
| 203 | + | |
| 204 | +\end{document} | |
| 205 | + | ... | ... |