Spectrum Analysis

The differential flux \phi(E) per area, time and energy interval is defined as \phi(E) = \frac{dN}{dA\cdot dt\cdot dE}

Often \phi(E) is also referred to as \frac{dN}{dE} as observation time and effective collection area is a constant. The effective area is then defined as A_\textrm{eff}(E)=\epsilon(E)\cdot A_0. Note that at large distances R_0 the efficiency \epsilon(R_0) vanishes, so that the effective area is an (energy dependent) constant while A_0=\pi R_0^2 and the efficiency \epsilon(E) are mutually dependent.

For an observation with an effective observation time \Delta T=\sum\delta t_i, this yields in a given Energy interval \Delta E: \phi(\Delta E) = \frac{1}{A_0\cdot \Delta T}\frac{N(\Delta E)}{\epsilon(\Delta E)\cdot \Delta E}

For simplicity, in the following, \Delta E will be replaced by just E but always refers to to a given energy interval. If data is binned in a histogram, the relation between the x-value for the bin E_\textrm{x} and the corresponding interval \Delta E_\textrm{x} is not well defined. Resonable definitions are the bin center (usually in logarithmic bins) or the average energy.

The total area A_0 and the corresponding efficiency \epsilon(E) are of course only available for simulated data. For simulated data, A_0 is the production area and \epsilon(E) the corresponding energy dependent efficiency of the analysis chain. For a given energy bin, the efficiency is then defined as

\epsilon(E) = \frac{N_\textrm{exc}(E)}{N_0(E)}

where N_0 is the number of simulated events in this energy bin and N=N_{exc} the number of *excess* events that are produced by the analysis chain.

Note that the exact calculation of the efficiency \epsilon(\Delta E) depends on prior knowledge of the correct source spectrum N_0(E). Therefore, it is strictly speaking only correct if the simulated spectrum and the real spectrum are identical. As the real spectrum is unknown, special care has to be taken of the systematic introduced by the assumption of N_0(E).

Zenith Angle Weights and Spectral Weights

While the source spectrum N_0 is of course independent of the zenith angle \theta of the observation, the efficiency is not, so that

\epsilon(E) = \epsilon(E, \theta)

In addition, the observed differential flux depends on the zenith angle. That means that for a correct calculation of the efficiency, the number of simulated events per zenith angle interval has to match the observation time distribution.

\epsilon(E, \theta) = \epsilon(E, \theta, \Delta T(\theta))

This can be achieved by applying zenith angle dependent weights \omega(\theta) to each simulated event. Similarly, the discussed requirement of the match between result spectrum and simulated spectrum can be achieved applying energy dependent weights \omega(E).

As the simulated energy spectrum is independent of zenith angle, it can be expressed as

\phi(E,\theta) = N_0\cdot \eta(E)\cdot \eta(\theta)

with the differential energy spectrum \eta(E) and the zenith angle distribution \eta(\theta) and the normalization

\int_E\int_\theta\phi(E,\theta)dEd\theta = N_0\cdot \int_E\eta(E)dE\cdot \int_\theta\eta(\theta)d\theta

The contents of one energy bin \Delta E and one zenith angle bin \Delta\theta can then be written as

N(\Delta E, \Delta\theta) = \sum_{\Delta E}\sum_{\Delta\theta} \omega_{i,j}(E, \theta)= \sum_{\Delta E}\sum_{\Delta\theta} \omega_i(E)\cdot\omega_j(\theta)= \sum_{\Delta E} \omega_i(\theta)\cdot\sum_{\Delta\theta}\omega_i(E)

and the normalization

\sum_E\sum_\theta N(E,\theta) = N_0 \sum_E\sum_\theta \omega_{i,j}(E,\theta) = N_0 \sum_E\sum_\theta \omega_{i}(E)\omega_j(\theta) = N_0\cdot \sum_E\omega_i(\theta)\cdot \sum_\theta\omega_i(E)

Excess

The number of excess events, for data and simulations, is defined as

N_\textrm{exc} = N_\textrm{sig} - \hat N_\textrm{bg}

where N_\textrm{sig} is the number of events identified as potential gammas from the source direction ('on-source') and N_\textrm{bg} the number of gamma-like events measured 'off-source'. Note that for Simulations, \hat N_\textrm{bg} is not necessarily zero for wobble-mode observations as an event can survive the analysis for on- and off-events, if this is not prevented by the analysis (cuts).

The average number of background events \hat N_\textrm{bg} is the total number of background events N_\textrm{bg} from all off-regions times the corresponding weight \frac{1}{5} (often referred to as \alpha). For five off-regions, this yields

\hat N_\textrm{bg} = \frac{N_\textrm{bg}}{5}

Assuming Gaussian errors, the statistical error is thus

\sigma^2(N_\textrm{exc}) = \left(\frac{dN_\textrm{exc}}{dN_\textrm{sig}}\right)^2\sigma^2(N_\textrm{sig}) + \left(\frac{dN_\textrm{exc}}{d\hat N_\textrm{bg}}\right)^2\sigma^2(\hat N_\textrm{bg})

For data this immediately resolves to

\sigma^2(N_\textrm{exc}) = N_\textrm{sig} + \frac{1}{5^2}N_\textrm{bg}

with the Poisson (counting) error \sigma^2(N_\textrm{sig,bg}) = N_\textrm{sig,bg}.

For the simulations, the error on N is the error on the weighted sum

\sigma^2(N_\textrm{sig,bg}) = \sum_i\left[\frac{\sum_j d\omega_j}{d\omega_i}\sigma(\omega_i)\right]^2 = \sum_i\sigma^2(\omega_i)

<!--With the error \sigma^2(\omega_i) = \omega_i^2 this provides \sigma^2(N_\textrm{exc}^\textrm{MC}) = \sum_\textrm{sig}\omega_i^2 + \frac{1}{5^2}\sum_\textrm{bg}\omega_i^2 -->

Code

In the following N refers to a number of simulated events and M to a number of measured (excess) events.

n(\delta E, \delta\theta) is the number of produced events in the energy interval E\in\delta E=[e_\textrm{min};e_\textrm{max}] and the zenith angle interval \theta\in\delta\theta=[\theta_\textrm{min};\theta_\textrm{max}]. The weighted number of events n'(\delta E, \delta\theta) in that interval is then

n'(\delta E, \delta\theta) = \sum_{i=0...n}^{E\in\delta E}\sum_{j=0...n}^{\theta\in\delta\theta} \omega(E_i, \theta_j)

Since \omega(E, \theta) = \rho(E)\cdot \tau(\theta) with \rho(E) the spectral weight to adapt the spectral shape of the simulated spectrum to the real (measured) spectrum of the source and \tau(\theta) the weight to adapt to oberservation time versus zenith angle.

n'(\delta E, \delta\theta) = \sum_{i=0...n}^{E\in\delta E}\sum_{j=0...n}^{\theta\in\delta\theta} \rho(E_i)\cdot\tau(\theta_j) = \sum_{i=0...n}^{E\in\delta E}\rho(E_i)\cdot\sum_{j=0...n}^{\theta\in\delta\theta} \tau(\theta_j)

The weighted number of produced events N'(\Delta E, \Delta\Theta) in the total energy interval E\in\Delta E=[E_\textrm{min};E_\textrm{max}] and the total zenith angle interval \theta\in\Delta\Theta=[\Theta_\textrm{min};\Theta_\textrm{max}] is then

N'(\Delta E, \Delta\Theta) = \sum_{i=0...N}^{E\in\Delta E}\sum_{j=0...N}^{\theta\in\Delta\Theta} \rho(E_i)\cdot\tau(\theta_j) = \sum_{i=0...N}^{E\in\Delta E}\rho(E_i)\cdot\sum_{j=0...N}^{\theta\in\Delta\Theta} \tau(\theta_j)

with N(\Delta E, \Delta\Theta) being the total number of produced events.

The weights are defined as follow:

\rho(E) = \rho_0\frac{\phi_\textrm{src}(E)}{\phi_0(E)}

where \phi_0(E) is the simulated spectrum and \phi_\textrm{src} the (unknown) real source spectrum. \rho_0 is a normalization constant. The zenith angle weights \tau(\delta\theta) in the interval \delta\theta are defined as

\tau(\delta\theta) = \tau_0\frac{\Delta T(\delta\theta)}{N(\delta\theta)}

where N(\delta\theta) is the number of produced events in the interval \delta\theta and \Delta T(\delta\theta) is the total observation time in the same zenith angle interval. \tau_0 is the normalization constant.

As the efficiency \epsilon(\delta E) for an energy interval \delta E is calculated as

\epsilon(\delta E) = \epsilon(\delta E, \Delta\Theta) = \frac{N'_\textrm{exc}(\delta E, \Delta\Theta)}{N'_\textrm{src}(\delta E,\Delta\Theta)}

and N'_\textrm{exc}(\delta E,\Delta\Theta) and N'_\textrm{src}(\delta E,\Delta\Theta) are both expressed as the sum given above, the constants \rho_0 and \tau_0 cancel.

The differential flux in an energy interval \delta E is then given as

\phi(\delta E) = \phi(\delta E,\Delta\Theta) = \frac{1}{A_0\cdot \Delta T}\frac{M_\textrm{exc}(\delta E)}{\epsilon(\delta E)\cdot \delta E} = \frac{M_\textrm{exc}(\delta E)}{N'_\textrm{exc}(\delta E)}\cdot \frac{N'_\textrm{src}(\delta E)}{A_0\cdot\Delta T\cdot \delta E}

Where A_0 is total area of production and \Delta T the total observation time. The number of measured excess events is in that energy interval is M_\textrm{exc}(\delta E)=M_\textrm{exc}(\delta E, \Delta\Theta).

Using Gaussian error propagation, the error in a given energy interval \delta is then given by

\sigma^2(\phi) = \left(\frac{1}{A_0\cdot \Delta T\cdot\delta E}\right)^2\cdot\left[\left(\frac{d\phi'}{dM_\textrm{exc}}\right)^2\sigma^2(M_\textrm{exc}) + \left(\frac{d\phi'}{dN'_\textrm{exc}}\right)^2\sigma^2(N'_\textrm{exc}) + \left(\frac{d\phi'}{dN'_\textrm{src}}\right)^2\sigma^2(N'_\textrm{src})\right]

with \phi':=A_0\cdot\Delta T\cdot\delta E\cdot\phi.

\frac{d\phi'}{dM_\textrm{exc}} = \frac{N'_\textrm{src}}{N'_\textrm{exc}} \frac{d\phi'}{dN'_\textrm{exc}} = -\frac{M_\textrm{exc}\cdot N'_\textrm{src}}{\left(N'_\textrm{exc}\right)^2} \frac{d\phi'}{dN'_\textrm{src}} = \frac{M_\textrm{exc}}{N'_\textrm{exc}}

\sigma^2(\phi) = \left(\frac{1}{A_0\cdot \Delta T\cdot\delta E}\right)^2\cdot\left[\left(\frac{N'_\textrm{src}}{N'_\textrm{exc}}\right)^2\sigma^2(M_\textrm{exc}) + \left(\frac{M_\textrm{exc}\cdot N'_\textrm{src}}{(N'_\textrm{exc})^2}\right)^2\sigma^2(N'_\textrm{exc}) + \left(\frac{M_\textrm{exc}}{N'_\textrm{exc}}\right)^2\sigma^2(N'_\textrm{src})\right]

= \left(\frac{1}{A_0\cdot \Delta T\cdot\delta E\cdot N'_\textrm{exc}}\right)^2\cdot\left[\left(N'_\textrm{src}\right)^2\sigma^2(M_\textrm{exc}) + \left(\phi'\right)^2\sigma^2(N'_\textrm{exc}) + \left(M_\textrm{exc}\right)^2\sigma^2(N'_\textrm{src})\right]

\sigma^2(N') = \sum_m\left(\frac{dN'}{d\rho_m}\sigma(\rho_m)\right)^2+\sum_n\left(\frac{dN'}{d\tau_n}\sigma(\tau_n)\right)^2

=\sum_m\left(\sum_j\tau(\theta_j)\sigma(\rho_m)\right)^2+\sum_n\left(\sum_i\rho(E_i)\sigma(\theta_n)\right)^2

=\sum_m\sigma^2(\rho_m)\left(\sum_j\tau(\theta_j)\right)^2+\sum_n\sigma^2(\theta_n)\left(\sum_i\rho(E_i)\right)^2

\sigma^2(\rho_m) = 0

\sigma^2(\theta_n) = \frac{\tau_0^2\Delta T_n^2}{N_n^4}\sigma^2(N_n)= \frac{\tau_0^2\Delta T_n^2}{N_n^3}

\sigma^2(N') = \tau_0^2\left(\sum_{n\in\Delta\Theta}\frac{\Delta T_n^2}{N_n^3}\right)\left(\sum_{i\in\delta E}\rho(E_i)\right)^2

Monte Carlo

\omega_i(E) = \frac{\phi_\textrm{src}(E)}{\phi_0(E)}=\frac{E^{-\gamma}}{E^{-2.7}}

\omega_i(\Delta \theta) = \frac{\sum_{\Delta\theta}\Delta T_i}{\sum_\theta\Delta T_i} \cdot \frac{M_0}{M_0(\Delta \theta)}

\epsilon(E) = \frac{N_\textrm{exc}(E)}{N_0(E)}

\phi(E) = \frac{N_\textrm{exc}(E)}{A_\epsilon(E)\cdot\Delta T} = \frac{N_\textrm{exc}(E)}{N_\textrm{exc}^\textrm{MC}(E)}\cdot \frac{N_0(E)}{A_0\cdot\Delta T}

\phi(E) = \frac{N_\textrm{sig}(E) - \hat N_\textrm{bg}(E)}{N_\textrm{sig}^\textrm{MC}(E) - \hat N_\textrm{bg}^\textrm{MC}(E)}\cdot \frac{N_0(E)}{A_0\cdot\Delta T}

N_\textrm{sig}^\textrm{MC} = \sum_\textrm{sig}^\textrm{MC}\omega_i(E)\omega_i(\theta) N_\textrm{bg}^\textrm{MC} = \sum_\textrm{bg}^\textrm{MC}\omega_i(E)\omega_i(\theta) N_0(E) = \sum_\textrm{corsika}^\textrm{MC}\omega_i(E)\omega_i(\theta)

\sigma^2(\phi(E)) = \left(\frac{d\phi(E)}{dN_0}\right)^2\sigma^2(N_0) + \left(\frac{d\phi(E)}{dN_\textrm{exc}^\textrm{MC}}\right)^2\sigma^2(N_\textrm{exc}^\textrm{MC}) + \left(\frac{d\phi(E)}{dN_\textrm{exc}}\right)^2\sigma(N_\textrm{exc})^2

\sigma^2(N_0) = \sum_\textrm{corsika}\omega_i^2(E)\omega_i^2(\theta)

\sigma^2(N_\textrm{exc}^\textrm{MC}) = \left(\frac{dN_\textrm{exc}^\textrm{MC}}{dN_\textrm{sig}^\textrm{MC}}\right)^2\sigma^2(N_\textrm{sig}^\textrm{MC}) + \left(\frac{dN_\textrm{exc}^\textrm{MC}}{d\hat N_\textrm{bg}^\textrm{MC}}\right)^2\sigma^2(\hat N_\textrm{bg}^\textrm{MC})

\sigma^2(N_\textrm{sig}) = N_\textrm{sig}

\sigma^2(\hat N_\textrm{bg}) = \frac{1}{5^2}\sigma^2(N_\textrm{bg}) = \frac{1}{5^2} N_\textrm{bg}

\sigma^2(N_\textrm{sig}^\textrm{MC}) = \sum_\textrm{sig}^\textrm{MC}\omega^2_i(E)\omega^2_i(\theta)

\sigma^2(\hat N_\textrm{bg}^\textrm{MC}) = \frac{1}{5^2}\sigma^2(N_\textrm{bg}^\textrm{MC}) = \frac{1}{5^2}\sum_\textrm{bg}^\textrm{MC}\omega^2_i(E)\omega^2_i(\theta)

\frac{dN_\textrm{exc}}{dN_\textrm{sig}}=1 \frac{dN_\textrm{exc}}{d\hat N_\textrm{bg}}=1

Errors

For a sum of weights, e.g. R = \sum_N \rho_i\tau_i the corresponding error is \[\sigma_R = \sum_N\left[\rho_i\sigma(\tau_i)+\tau_i\sigma(\rho_i)\right]²

Theta and Spectral Weights

While the source spectrum N_0 is of course independent of the zenith angle \theta of the observation, the efficiency is not, so that

\epsilon(E) = \epsilon(E, \theta)

In addition, the observed differential flux depends on the zenith angle, so that

\epsilon(E, \theta) = \epsilon(E, \theta, \Delta T(\theta))

That means that for a correct calculation of the efficiency, the number of simulated events per zenith angle interval has to match the observation time distribution. This can be achieved by applying zenith angle dependent weights \omega(\theta) to each simulated event. Similarly, the discussed requirement of the match between result spectrum and simulated spectrum can be achieved applying energy dependent weights \omega(E).

As the simulated energy spectrum is independent of zenith angle, it can be expressed as

N_0(E,\theta) = N_0\cdot \fract{\eta(E)}{\int_E\eta(E)dE}\cdot \fract{\eta(\theta)}{\int_\theta\eta(\theta)d\theta}

with the differential energy spectrum \eta(E) and the zenith angle distribution \eta(\theta).

dN_0(E,\theta) = N_0\cdot d\eta(E)\cdot d\eta(\theta)

N_0 = \int_{E_\textrm{min}}^{E_\textrm{max}}\int_{\theta_\textrm{min}}^{\theta_\textrm{max}} dN_0(E,\theta) = N_0\int_{E_\textrm{min}}^{E_\textrm{max}}\int_{\theta_\textrm{min}}^{\theta_\textrm{max}}\eta(E) dE \cdot \eta(\theta) d\theta= N_0\int_{E_\textrm{min}}^{E_\textrm{max}}\eta(E) dE\cdot \int_{\theta_\textrm{min}}^{\theta_\textrm{max}}\eta(\theta) d\theta

consequently

\int_{E_\textrm{min}}^{E_\textrm{max}}\eta(E) dE = \int_{\theta_\textrm{min}}^{\theta_\textrm{max}}\eta(\theta) d\theta = 1

or

N_0 = N_0\cdot \sum_{E_\textrm{min}}^{E_\textrm{max}}\omega_i(E) \cdot \sum_{\theta_\textrm{min}}^{\theta_\textrm{max}}\omega_i(\theta)

with

\sum_{E_\textrm{min}}^{E_\textrm{max}}\omega_i(E) =\sum_{\theta_\textrm{min}}^{\theta_\textrm{max}}\omega_i(\theta) = 1

\int_{E_\textrm{min}}^{E_\textrm{max}}\int_{\theta_\textrm{min}}^{\theta_\textrm{max}}\eta(E)\eta(\theta) dE d\theta =\int_{E_\textrm{min}}^{E_\textrm{max}}\eta(E) dE = \int_{\theta_\textrm{min}}^{\theta_\textrm{max}}\eta(\theta) d\theta = 1

and N_0 the total number of generated Monte Carlo events. For the generated number of events n_0 in the energy interval \Delta E=E_\textrm{max}-E_\textrm{min} and zenith angle interval \Delta \theta=\theta_\textrm{max}-\theta_\textrm{min} this is

n_0(\Delta E) = \frac{\sum_{\Delta E}\sum_{\theta_\textrm{min}}^{\theta_\textrm{max}} \omega_i(E)\cdot\omega_i(\theta)}{\sum_{E}\sum_{\theta} \omega_i(E)\omega_i(\theta)}

n_0(\Delta E) = \frac{\sum_{\Delta E}\sum_{\theta_\textrm{min}}^{\theta_\textrm{max}} \omega_i(E)\cdot\omega_i(\theta)}{\sum_{\Delta E}\sum_{E_\textrm{min}}^{E_\textrm{max}} \omega_i(E) \cdot \sum_{\Delta \theta}\sum_{\theta_\textrm{min}}^{\theta_\textrm{max}} \omega_i(\theta)}

with the weights chosen such that the sum over all intervals

\sum_{\Delta E}\sum_{E_\textrm{min}}^{E_\textrm{max}} \omega_i(E)=\sum_{\Delta \theta}\sum_{\theta_\textrm{min}}^{\theta_\textrm{max}} \omega_i(\theta)=1

Define Binnings

Get Data File List

Get Observation Time

Get Monte Carlo File List

Get Zenith Angle Histogram

Analyze Data

Analyze Monte Carlo Data

Summarize Corsika Production

Result (Spectrum)

Result (Threshold)

Result (Migration)