= 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 the efficiency vanishes, so that the effective area is an (energy dependent) constant while \(A_0\) and the efficiency \(\epsilon(E)\) are mutually dependent.

For an observation with an effective observation time \(\Delta T\), this yields in a given Energy interval \(\Delta E\):
\[\phi(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. 

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\). 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\).

=== 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 protected by the analysis chain.

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 \(\omega\). For five off-regions, this yields

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

=== Monte Carlo ===

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

\[\phi(E) = \frac{N_\textrm{exc}(E)}{A_0\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)\sigma^2(N_0) + \left(\frac{d\phi(E)}{dN_\textrm{exc}^\textrm{MC}}\right)\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{d\phi(E)}{dN_\textrm{sig}^\textrm{MC}}\right)\sigma^2(N_\textrm{sig}^\textrm{MC}) +  \left(\frac{d\phi(E)}{d\hat N_\textrm{bg}^\textrm{MC}}\right)\sigma^2(\hat N_\textrm{bg}^\textrm{MC}) \]
\[\sigma^2(N_\textrm{exc}) = \left(\frac{d\phi(E)}{dN_\textrm{sig}}\right)\sigma^2(N_\textrm{sig}) +  \left(\frac{d\phi(E)}{d\hat N_\textrm{bg}}\right)\sigma^2(\hat N_\textrm{bg}) \]

\[\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}) = \]
\[\sigma^2(\hat N_\textrm{bg}^\textrm{MC}) = \frac{1}{5^2}\sigma^2(N_\textrm{bg}^\textrm{MC}) = \]





=== 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

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

with the differential energy spectrum \(\eta(E)\) and the zenith angle distribution \(\eta(\theta)\), for the total production range \(E_\textrm{min}\)  to \(E_\textrm{max}\) and \(\theta_\textrm{min}\) to \(\theta_\textrm{max}\)



\[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) ==