Awakening the BALROG: BAyesian Location Reconstruction Of GRBs

Awakening the BALROG: BAyesian Location Reconstruction Of GRBs Abstract The accurate spatial location of gamma-ray bursts (GRBs) is crucial for both accurately characterizing their spectra and follow-up observations by other instruments. The Fermi Gamma-ray Burst Monitor (GBM) has the largest field of view for detecting GRBs as it views the entire unocculted sky, but as a non-imaging instrument it relies on the relative count rates observed in each of its 14 detectors to localize transients. Improving its ability to accurately locate GRBs and other transients is vital to the paradigm of multimessenger astronomy, including the electromagnetic follow-up of gravitational wave signals. Here we present the BAyesian Location Reconstruction Of GRBs (balrog) method for localizing and characterizing GBM transients. Our approach eliminates the systematics of previous approaches by simultaneously fitting for the location and spectrum of a source. It also correctly incorporates the uncertainties in the location of a transient into the spectral parameters and produces reliable positional uncertainties for both well-localized sources and those for which the GBM data cannot effectively constrain the position. While computationally expensive, balrog can be implemented to enable quick follow-up of all GBM transient signals. Also, we identify possible response problems that require attention and caution when using standard, public GBM detector response matrices. Finally, we examine the effects of including the uncertainty in location on the spectral parameters of GRB 080916C. We find that spectral parameters change and no extra components are required when these effects are included in contrast to when we use a fixed location. This finding has the potential to alter both the GRB spectral catalogues and the reported spectral composition of some well-known GRBs. methods: data analysis, methods: statistical, gamma-ray burst: general 1 INTRODUCTION The localization of gamma-ray bursts (GRBs) on the sky has played an important role in the history of understanding these extreme objects. The Burst and Transient Source Experiment (BATSE) helped to confirm the cosmological origin of GRBs by determining the isotropic sky distribution of over 3000 GRBs (Briggs et al. 1996). The Fermi Gamma-ray Burst Monitor (GBM), which can monitor the entire sky (with the exception of the portion occulted by the Earth), is currently providing locations for hundreds of GRBs and other transient gamma-ray events every year. With the advent of gravitational wave astronomy and the potential for electromagnetic follow-up of associated events, GBM could play a crucial role in the new era of multimessenger astronomy. Using differential rates in each of its 14 non-imaging detectors, GBM is able to provide approximate localization of transient events that can then be used by imaging observatories to search for counterparts at other wavelengths. Unfortunately, the standard GBM localization procedure, the Daughter Of Location (DOL), produces offsets of 8°–13° from the subsequently confirmed true locations for many GRBs localized by GBM (Connaughton et al. 2015). The exact cause of this systematic is unknown and could be attributed to several factors such as imperfections in the simulated spacecraft response and/or model. Here we hypothesize that the systematics in the standard GBM localization algorithm are a by-product of the assumption that all GRB spectral energy distributions (SEDs) are alike, and of inaccuracies in the detector response matrix1 (DRM) modelling. The DOL uses a set of three template spectra are assumed in order to calculate a spatial grid of putative observed counts on the sky. This grid is then compared with the real observed counts via a χ2-minimization on source position to produce a best-fitting location. The template spectra are built from the ubiquitous Band function (Band, Matteson & Ford 1993) with three sets of fixed spectral parameters that allow for only the amplitude to be varied. Typically, one template is used to locate a GRB, while the other two are reserved for solar flares and soft-gamma repeaters. Since GRBs exhibit a variety of spectral shapes and temporal spectral evolution (e.g. Burgess 2014), the assumption of a fixed set of spectral shapes could influence inferred location of a GRB if its true SED differs from the templates. To explore this possibility we have developed a method where the location and spectrum are inferred simultaneously via Bayesian sampling using a new method: the BAyesian Location Reconstruction Of GRBs (balrog).2 First, we detail our methodology (Section 2) then demonstrate the viability of the method via simulations (Section 3) and real data (Section 4). Finally, we investigate the effect on spectra introduced by including the uncertainty of location where we use GRB 080916C as a case study (Section 5). 2 INFERENCE METHODOLOGY We analyse GBM data by using Bayesian inference to obtain the posterior distribution of the angular position and emission spectrum of a GRB. This requires a model of the GBM instrument (Section 2.1), a GRB source (Section 2.2) and the sky background (Section 2.3), which can then be combined to obtain the full likelihood for the data set (Section 2.4). Combining the likelihood with the parameter priors (Section 2.5), the implied constraints are then obtained by generating samples from the posterior distribution (Section 2.6). This allows for easy marginalization over nuisance parameters, most obviously integrating out the parameters describing the GRB's spectrum to focus on localization. 2.1 GBM instrumental response The on-ground calibration of GBM was described in Bissaldi et al. (2009), which was used as input for designing the GBM DRMs. Any method of localizing and characterizing sources using GBM must account for the complicated nature of its detectors, which have a broad angular response and suffer from energy dispersion (i.e. there is a significant probability that only a fraction of a detected photon's energy is registered), all of which should be embodied in the DRMs. These effects are accounted for using forward folding, described in detail by Vianello et al. (2015), which provides a probabilistic connection between the actual properties of incident photons and the data recorded by the GBM detectors. A critical simplifying approximation that helps make this problem tractable is that the orientations of the Nd = 14 GBM detectors are perfectly known and constant for the duration of the burst being analysed. This means that the rotation angles (represented as quaternions) that define each detector's orientation can be left implicit, although the formalism presented below could be extended to include the changing orientation of the detectors with the satellite's motion. During the on-source period of a GRB each detector is exposed to a radiation field that here is described by the number flux of photons at the detector (i.e. the rate per unit area, per unit energy, per steradian), n(E, ψ), which depends on both photon energy, E, and angular position, ψ, but is assumed to be time independent. The rate at which photons are registered by detector d (with d ∈ {1, 2, …, Nd}) is n(E, ψ) Ad(E, ψ), where Ad(E, ψ) is the effective area of the detector. The angular dependence of the effective area can largely be understood in terms of the geometrical area presented to a putative source, and so it decreases approximately with the cosine of the angle between ψ and nominal detector orientation, although the output of a full numerical simulation is used in practice. Each photon that is registered then deposits an energy ε ≲ E, a process which can be fully described by the conditional distribution Pr(ε|E, ψ, d, detected). The most important aspect of this energy dispersion is that there is a significant probability that ε is considerably lower than E, due to photons Compton scattering then escaping before all their energy is deposited in the scintillation crystal. This aspect of the detector response has an angular dependence, inextricably linking the probability that a given photon is detected with the energy recorded from it. In the GBM analysis formalism photon detection and energy registration are combined into a single function,   \begin{equation} R_d(\varepsilon , E, \psi ) = A_d(E, \psi ) \, {\rm Pr}(\varepsilon | E, \psi , d, {\rm detected}), \end{equation} (1)which is stored as a DRM and fully characterizes each detector. The energy of an individual photon is not recorded per se, but rather each photon is recorded as being in one of Nc = 128 channels that have a nominal correspondence with energy. The rate at which photons are recorded in channel c (with c ∈ {1, 2, …, Nc} ) is   \begin{equation} S_{d,c} = \int _{E_{{\rm min}, {\rm c}}}^{E_{{\rm max}, {\rm c}}} {\rm d}\varepsilon \int _{0}^\infty {\rm d}E \int {\rm d}\Omega \, n(E, \psi ) \, R_d(\varepsilon , E, \psi ), \end{equation} (2)where Emin, c and Emax, c are the minimum and maximum energies for this channel and the inner integral extends over all photon arrival directions, ψ. 2.1.1 DRM generation GBM DRMs are built via the General Response Simulation System (gress) geant4 Monte Carlo code (Kippen et al. 2007) and separated into two components consisting of a direct DRM and an atmospheric scattering DRM of the 14 detectors for each individual GRB position. A grid of direct DRMs created for 272 points on the sky for each of the 12 GBM sodium iodide (NaI) detectors and two bismuth germinate (BGO) detectors is used as the base response library for DRM generation. These responses exist in an energy compressed format that must be decompressed and added to the atmospheric and satellite internal scattering response via a series of transformations that depend on GBM's position and orientation with respect to the Earth (Paciesas et al. 1988). The directional dependence of the response adds in an extra complication, as it is too computationally expensive to calculate Rd(ε, E, ψ) for all values of ψ over the sphere, so this operation is undertaken only for the source positions proposed during the posterior sampling process (Section 2.6). Therefore, we designed efficient DRM generation software based on the original GBM tool that allows for vectorized and parallel generation of DRMs. The final direct DRM is added to the atmospheric scattering DRM to account for scattering of photons off the Earth's atmosphere (Meegan et al. 2009). balrog has the ability to correctly reject source positions that are occulted by the Earth even though these positions are accounted for in the base DRM data base. The area of the sky that is occulted is calculated in real time using the current orbital altitude of the Fermi satellite and finding the horizon line from the spacecraft zenith. But it is more computationally efficient to do this before generating the full detector response, so the DRM is immediately set to Rd(ε, E, ψ) = 0 for occulted directions. 2.2 Source model For an astronomically distant point source (as GRBs are assumed to be) all photons arrive at the detector from the direction of the source, ψs. This allows the number flux defined in Section 2.1 to be factorized as  \begin{equation} n(E, \psi ) = f(E, \,\,{\phi }_{\rm s}) \, \delta _{\rm D}(\psi , \psi _{\rm s}), \end{equation} (3)where $$\,\,{\phi }_{\rm s}$$ is the vector of parameters (normalization, spectral index, etc.) describing the source's SED and δD(ψ, ψs) is the delta function distribution on the sphere centred on ψs. The presence of the delta function means the predicted count rate from the source given in equation (2) simplifies to be   \begin{equation} S_{d,c} (\psi _{\rm s}, \,\,{\phi }_{\rm s}) = \int _{0}^\infty {\rm d}E \, f(E, \,\,{\phi }_{\rm s}) \int _{E_{{\rm min}, {\rm c}}}^{E_{{\rm max}, {\rm c}}} {\rm d}\varepsilon \, R_d(\varepsilon , E, \psi _{\rm s}), \end{equation} (4)where the dependence on the source parameters has been made explicit. The order of integration has been swapped to emphasize the possibility of pre-calculating the inner integral independent of the source emission parameters. The SED models that are used here are Band, cut-off power law (CPL) and smoothly broken power law (SBPL). The SED parameters $$\,\,{\phi }_{\rm s}$$ and the source position ψs represent the full set of GRB parameters to be constrained using the GBM data. 2.3 Background estimation GBM data include an isotropic sky background and an anisotropic Earth albedo background, both of which are time- and energy dependent. Accounting for the combined background rate from these two effects is a necessary step to inferring the properties of the GRBs detected by the instrument. The ideal approach to this problem would be to simultaneously model the background and source (cf. Loredo & Lamb 1992; van Dyk, Connors & Kashyap 2001), although the fact that the background counts are generally fairly high means that it is a good approximation to estimate the background rate separately. The background is taken to be a polynomial in time, and the coefficients are fit to the off-source light curve in each energy channel, with a Poisson likelihood used to ensure valid results even for high-energy channels with low counts. The result of this process is a full joint posterior distribution of the polynomial coefficients. This posterior is then propagated to the on-source period, the result of which is a data-driven estimate of the background rate, $$\hat{B}_{d,c}$$, and its uncertainty, σd, c, in channel c of detector d. Only the former would be used if a simple background subtraction were sufficient; but, as has been shown by Greiner et al. (2016), failure to include the background uncertainty can lead to qualitatively incorrect inferences about the source of interest. The implied posterior distribution for the background is taken to be a truncated normal of the form   \begin{eqnarray} {\rm Pr}(B_{d,c} | \hat{B}_{d,c}, \sigma _{d,c}) & = & \frac{\Theta (B_{d,c}) }{\left[1 + {\rm erf}\left(\hat{B}_{d,c} / 2^{1/2} \, \sigma _{d,c}\right) \right] / 2}\\ & &\times \frac{1}{(2 \pi )^{1/2} \, \sigma _{d,c}}\, {}{\rm e}^{-\left(B_{d,c} - \hat{B}_{d,c}\right)^2 / \left(2 \, \sigma ^2_{d,c} \right)},\nonumber \end{eqnarray} (5)where erf( · ) is the standard error function and the Heaviside step function Θ( · ) imposes the absolute prior constraint that the background rate cannot be negative. 2.4 Likelihood The likelihood, $${\rm Pr}(\lbrace N_{d,c}\rbrace | \psi _{\rm s}, \,\,{\phi }_{\rm s})$$, is the probability, conditional on a given source position ψs and emission model $$\,\,{\phi }_{\rm s}$$, of obtaining the measured counts, {Nd, c}, during the period (of duration T) that the GBM was observing the GRB. As such, it encodes all the relevant information that can be used to localize and characterize the source. Assuming that each incident photon is only registered once (i.e. in a particular channel of a particular detector), the likelihood has the form   \begin{equation} {\rm Pr}(\lbrace N_{d,c}\rbrace | \psi _{\rm s}, \,\,{\phi }_{\rm s}) = \prod _{d = 1}^{N_{\rm d}} \prod _{c = 1}^{N_{\rm c}} {\rm Pr}\left[N_{d,c} | S_{d,c} (\psi _{\rm s}, \,\,{\phi }_{\rm s})\right], \end{equation} (6)where $$S_{d,c} (\psi _{\rm s}, \,\,{\phi }_{\rm s})$$ is the predicted rate at which source photons are registered in channel c of detector d is given in equation (4). 2.4.1 Single channel contribution As the contribution from each channel and detector combination has the same mathematical form, it is simplest to consider this detector/channel combination without reference to detector number d or channel index c, which are hence left implicit in what follows. If the background rate, B, were known precisely, then the probability of N counts being registered during the on-source period would be given by a simple Poisson distribution as   \begin{equation} {\rm Pr}(N | S, B) = \frac{[(S + B)\, T]^{N} \, {\rm e}^{- (S + B)\, T}}{N!}, \end{equation} (7)where S is the rate of photons being registered from the source. This depends on ψs and $$\,\,{\phi }_{\rm s}$$, but it can be considered as the one model parameter of interest here. The fact that the background is not known perfectly must be incorporated, which formally should be done by including B as an extra parameter to be subsequently marginalized out. Fortunately, however, almost all the available information about B comes from the polynomial fitting to the off-source intervals described in Section 2.3, and not the single on-source measurement. It is hence possible to simply average the Poisson distribution given in equation (7) over the plausible values of the background rate. Given the background estimate $$\hat{B}$$ and associated uncertainty σ, and adopting the truncated normal defined in equation (5), the contribution to the likelihood from a single channel is modified to be   \begin{eqnarray} {\rm Pr}(N | S, \hat{B}, \sigma ) & = & \int _0^\infty {\rm d}B \, {\rm Pr}(B | \hat{B}, \sigma )\, {\rm Pr}(N | S, B)\nonumber \\ & \propto & \int _0^\infty {\rm d}B \, (S + B)^{N} {}{\rm e}^{- [(B + S)\, T + (B - \hat{B})^2 / (2 \, \sigma ^2)]},\nonumber \\ \end{eqnarray} (8)where only the terms with that have a contribution from the source parameters have been retained in the second expression. While the above integral could be evaluated numerically, the fact that it will be calculated numerous times during the posterior sampling process motivates using a closed form approximation instead. The option taken here is to evaluate the profile likelihood, in which, for a given value of S, the integral is taken to be the value of the integrand evaluated at its peak (i.e. for the value of the background rate, $$\tilde{B}$$, that maximizes the integrand). This is the approach implemented by Arnaud, Dorman & Gordon (2015) to obtain the pgstat statistic that is an option in XSPEC. This is equivalent3 to assuming that   \begin{eqnarray} {\rm Pr}(N | S, \hat{B}, \sigma ) \propto \left\lbrace \begin{array}{ll}{\rm e}^{-[ (S + \hat{B})\, T - \sigma ^2 \, T^2 / 2 ]}, & \!\!\!\! {\rm if}\ N = 0, \\ (S + \tilde{B})^{N} \, {}{\rm e}^{- [(S + \tilde{B}) \, T + (\tilde{B} - \hat{B})^2 / (2 \, \sigma ^2)]}, & \!\! {\rm if}\ N > 0, \end{array}\right. \nonumber\\ \end{eqnarray} (9)where in the N > 0 case   \begin{eqnarray} \tilde{B} &=& \frac{1}{2} \! \left\lbrace \hat{B} -S -\sigma ^2 T + \! \left[ (S + \sigma ^2 T - \hat{B})^2\right.\right.\nonumber\\ &&\left.\left. -\, 4 (\sigma ^2 T S - \sigma ^2 N - \hat{B} S) \right]^{1/2} \right\rbrace \! . \end{eqnarray} (10) Independent of whether numerical integration or the profile likelihood is used to evaluate $${\rm Pr}(N | S, \hat{B}, \sigma )$$, there is an increase in the range of values of S for which the likelihood is appreciable. This can lead to the apparently undesirable result that the parameter constraints are similarly broadened, but in fact this is simply a correct reflection of the range of source positions that are consistent with the data and hence which should be considered as potential GRB locations. 2.5 Parameter priors A necessary component of Bayesian parameter estimation is the specification of priors for all parameters. The joint prior in ψs and $$\,\,{\phi }_{\rm s}$$ can reasonably be assumed to factorize (i.e. the plausible SED parameters are independent of position on the sky), so that $${\rm Pr}(\psi _{\rm s}, \,\,{\phi }_{\rm s} | {\rm GRB}) = {\rm Pr}(\psi _{\rm s} | {\rm GRB}) \, {\rm Pr}(\,\,{\phi }_{\rm s} | {\rm GRB})$$. For ψs we assume a prior that is uniform over the celestial sphere, so that Pr(ψs|GRB) = 1/(4π), although there is no operational need to normalize this distribution correctly. For $$\,\,{\phi }_{\rm s}$$ we assume uniform improper priors with appropriate bounds determined for the spectral model at hand (e.g. to ensure that the SED is not negative for positive E). The choice of prior can of course be more physically motivated, however, for our purposes these choices are sufficient because, as we show in Section 4, the marginalized positional posteriors do not have significant prior dependence. 2.6 Posterior distribution and sampling Given the likelihood described in Section 2.4 and the parameter priors given in Section 2.5, the desired posterior distribution has the form   \begin{equation} {\rm Pr}(\psi _{\rm s}, \,\,{\phi }_{\rm s} | \lbrace N_{d,c}\rbrace , {\rm GRB}) \propto {\rm Pr}(\psi _{\rm s}, \,\,{\phi }_{\rm s} | {\rm GRB}) \, {\rm Pr}(\lbrace N_{d,c}\rbrace | \psi _{\rm s}, \,\,{\phi }_{\rm s}).\!\!\! \end{equation} (11)There is no analytic way of determining the combinations of parameter values for which the posterior is high, nor of obtaining the normalizing constant that is needed to marginalize over $$\,\,{\phi }_{\rm s}$$. As the posterior can be multimodal, we use the form of nested sampling (Skilling 2004) as implemented in multinest (Feroz, Hobson & Bridges 2009) to obtain samples from the posterior. A particular utility of representing the posterior distribution as a set of samples is that marginalization is straightforward. In our case, depending on the situation, either $$\,\,{\phi }_{\rm s}$$ or ψs are nuisance parameters and marginalizing over one or other allows us to examine the marginal distribution of the parameter(s) of interest, Pr(ψs|{Nd, c}, GRB) or $${\rm Pr}(\,\,{\phi }_{\rm s} | \lbrace N_{d,c}\rbrace , {\rm GRB})$$, respectively. In order to be useful for timely follow-up, the balrog must be able to localize a transient swiftly. The current implementation can be run on GBM data types including the quick-look TRIGDAT that consists of light curves for each of the 14 GBM detectors each with the nominal 128-channel data rebinned to 8-channel energy resolution, as well as GBM time-tagged event (TTE) data that consist of photon time tags with 128-channel energy resolution (Meegan et al. 2009). Localization on a multicore desktop workstation takes ∼5–10 min for TRIGDAT data and ∼10–20 min for TTE data due to its finer energy resolution. Faster speeds can be achieved with a computing cluster-based installation that is used in this study. We discuss the differences between results with TRIGDAT and TTE data in the following sections. 3 SYNTHETIC LOCATION ANALYSIS To assess the validity of our method, we generate synthetic GRBs with known location and spectrum on a 25-point grid on the sky for a given spacecraft orientation. The CPL SED is used for the synthetic photon spectrum with power-law index of −1 and νFν peak of 300 keV. For each point in the grid, a DRM is generated for the given sky position and the SED is convolved with the DRM and given Poisson variability to produce an 8-channel synthetic data spectrum resembling TRIGDAT data. A simple Poisson background is added to the data obeying a power-law energy spectrum. Using balrog, we fit each simulated GRB in the grid with the CPL model to determine the location and spectrum using all 12 NaI detectors. Fig. 1 shows the distribution of angular separations between the true simulated position and the best-fitting > position. This is always less than 1°and has a mean of ∼0$$_{.}^{\circ}$$1. For real sources, the error could be much larger and depends on the intensity of the GRB. These results depend on knowing the true model of the spectrum and the calibration of the NaI detectors being perfect. As we will see in the following sections, errors in the response model will lead to systematic deviations from the true location when all 12 NaI detectors are used. Figure 1. View largeDownload slide Distribution of angular separations between the true simulated positions and the best-fitting ? position using the CPL SED for the 25-point simulated grid. Figure 1. View largeDownload slide Distribution of angular separations between the true simulated positions and the best-fitting ? position using the CPL SED for the 25-point simulated grid. In order to demonstrate the effect of using fixed templates to determine locations, we fit the simulated grid with a Band function with spectral parameters Ep = 1 MeV, α = 0 and β = −1.5 (i.e. the hard template). The amplitude and location are the only parameters allowed to vary in the fit. Fig. 2 shows that using the hard template introduces a systematic offset in the location of ∼13°, which is comparable to that identified by Connaughton et al. (2015). Moreover, we can examine the angular separation between the best-fitting positions and the location posterior (essentially the error region of the location) for one point on the grid and compare it with the angular separation between the true position and the location posterior. Fig. 3 shows that, when fitting with the CPL model, these two distributions are indistinguishable, which indicates that the error region of the fit is encompassing the true uncertainty in the data. However, while the error region obtained from the hard template is very precise, it is not accurate (i.e. does not encompass the true error region). Figure 2. View largeDownload slide Distribution of angular separations from the simulated positions to the best-fitting 4 position using the hard template SED for the 25-point simulated grid. Figure 2. View largeDownload slide Distribution of angular separations from the simulated positions to the best-fitting 4 position using the hard template SED for the 25-point simulated grid. Figure 3. View largeDownload slide Comparison between distributions of the angular separation of the location posterior and the best-fitting location (blue) and true location (green) for the CPL and hard template SEDs. Figure 3. View largeDownload slide Comparison between distributions of the angular separation of the location posterior and the best-fitting location (blue) and true location (green) for the CPL and hard template SEDs. We can look at how the differences in the locations found by the CPL and hard template models affect the DRM used for spectral analysis. As a toy example, we propagate the location posterior through the DRM generator to create a probabilistic DRM for both models. Fig. 4 illustrates the difference between these two DRMs. The main differences are at high and low energy, which can cause different slopes and hence different spectral parameters to be recovered. However, in the balrog scheme, the DRM is no longer a fixed quantity, but rather a probabilistic quantity weighted by the uncertainty in location. Figure 4. View largeDownload slide A probabilistic difference between the DRMs generated from the location posterior suing the hard template and the CPL SED. The colour scale indicates the difference in effective area in arbitrary units. This serves to illustrate the differences that can arise in the normal public DRMs when an improper model is used to generate DRMs. Figure 4. View largeDownload slide A probabilistic difference between the DRMs generated from the location posterior suing the hard template and the CPL SED. The colour scale indicates the difference in effective area in arbitrary units. This serves to illustrate the differences that can arise in the normal public DRMs when an improper model is used to generate DRMs. 3.1 Posterior checking To assess whether balrog provides GRBs position estimates of GRBs with realistic uncertainties, we need to perform model checking. The standard Bayesian approach is a posterior predictive check (PPC; see e.g. Gelman et al. 2014), in which simulated data sets for the GRB, drawn from the posterior distribution, are compared to the actual data on that source. Here, however, we are in the unusual situation of knowing the true locations of the GRBs from either because there accurate positions from other instruments such as Swift or because the data were simulated with known parameters. That allows us to bypass the data-simulation step and simply assess whether our posterior distributions are consistent with the GRBs’ actual positions. For a single GRB, with known sky position ψtrue, the key statistical quantity is the posterior density at this location, ρtrue = Pr(ψs = ψtrue|{Nd, c}). In cases where the posterior is singly peaked and symmetric this reduces to a measure based on the angular separation between the peak of the posterior and the true position, but is also valid if the posterior has significant degeneracies or is multimodal. The numerical value of ρtrue is not meaningful in isolation, as it is a probability density that is determined by e.g. the scale of the positional uncertainties. It must instead be evaluated in relation to the distribution of ρ values that would be obtained by drawing random positions from the posterior. A sample from this distribution can be obtained by drawing a position ψsim from the posterior distribution Pr(ψs|{Nd, c}) and then evaluating ρsim = Pr(ψs = ψsim|{Nd, c}) (i.e. the density at this simulated position). Repeating this process yields a set of samples from the distribution of ρ values from which ρtrue ought to have been drawn if the model were correct. If the model is badly wrong, then the true position will lie outside the region identified by the posterior and ρtrue will be lower than (almost all of) the simulated values ρsim. Conversely, if the uncertainties are overestimated, then ρtrue will be higher than the simulated values. This comparison can be made on a source-by-source basis, but for a set of NGRB GRBs the ρ values can be multiplied to obtain   \begin{equation} \rho _{\rm true} \equiv \prod _{i = 1}^{N_{\rm GRB}} \rho _{{\rm true},i}, \end{equation} (12)where ρi ≡ Pr(ψs = ψtrue, i|{Nd, c}, i) is evaluated from the posterior distribution for the ith GRB. Similarly, with ρsim, i, j being the jth simulated ρ value drawn from the posterior of the ith GRB, a sample from the expected distribution is given by the product   \begin{equation} \rho _{{\rm sim},j} \equiv \prod _{i = 1}^{N_{\rm GRB}} \rho _{{\rm sim},i,j}. \end{equation} (13) Applying this framework to the simulations, we obtain Fig. 5. As expected from simulations assuming an ideal detector response, the model is accurate at predicting future predictions (ρtrue is not too far to the left) without overfitting the data with the increased uncertainty of the balrog model (ρtrue is not too far to the right). Figure 5. View largeDownload slide The model checking distribution for the simulated GRB sets using the full approach of the balrog. Figure 5. View largeDownload slide The model checking distribution for the simulated GRB sets using the full approach of the balrog. 4 ANALYSIS OF REAL GRBS WITH KNOWN LOCATIONS Accessing the validity of our method depends strongly on our ability to precisely and accurately locate real GRBs that have been located by imaging instruments. We use the table of known GRB locations provided in Connaughton et al. (2015) that additionally allows us to compare our method to GRBs localized by the standard GBM method. We divide the analysis into two sections: GRBs well localized by GBM DOL and GRBs localized with a large systematic deviation by GBM DOL. First, we note a few details about fitting real data. In the above simulations, we assumed the DRMs are perfect and expect to recover the true location and spectral parameters. If there are calibration/simulation effects not properly modelled by the DRMs, then analysis of real data can produce suboptimal results leading to offsets in location and spectral parameters. Similarly, if the spectral model we choose to fit the GRBs does not correspond to the true spectral model of the source, we can expect the location to be affected. In real-time trigger situations, the user of balrog will not know the location a priori. As we will demonstrate in the following sections, the user must select a subset of detectors to obtain correct locations of sources. The GBM team has assembled a table of so-called legal detector sets to aid in choosing proper subsets of detectors for analysis (see Fig. 6; Meegan, team internal information). It will often be required to perform several location runs, with different spectral models and detector subsets to find the best location (see Section 7). balrog provides the marginal likelihood or evidence that can be used to compare between these iterations and choose the best location. It is currently unlikely that an automated system can perform these actions due to the complicated interaction between background, source interval selection, model choice and detector selection. Figure 6. View largeDownload slide Chart of GBM NaI detectors that see the same source based on a simulation of 10  000 GRBs made by the Fermi GBM team. Colour and number correspond to the number of coincident detections between a pair of detectors. It is intended to aid in choosing appropriate sets of detectors for spectral and location analysis. Figure 6. View largeDownload slide Chart of GBM NaI detectors that see the same source based on a simulation of 10  000 GRBs made by the Fermi GBM team. Colour and number correspond to the number of coincident detections between a pair of detectors. It is intended to aid in choosing appropriate sets of detectors for spectral and location analysis. 4.1 GRBs localized with GBM DOL resulting in a systematic deviation We examine 10 GRBs with systematic deviation in the DOL location from the true location to test if balrog can improve upon the result or provide appropriate location contours. We define a systematic deviation such that the angular separation between the DOL location and the true location exceeds the reported DOL 1σ error. We choose the 10 worst cases from Connaughton et al. (2015) that are listed in Table 1. For each of these GRBs, we try three spectral models (Band, CPL and SBPL) and two detector selections: all detectors and a subset of detectors. Table 1. Systematic location GRB sample. GRB  DOL RA  DOL Dec.  DOL Err  True RA  True Dec.  Separation angle  080725435  123.1  −23.1  2.2  121.7  −14.0  9.19  081101491  80.8  11.6  10.1  95.8  −0.1  18.94  090927422  67.6  −67.6  12.1  344.0  −71.0  27.37  100116897  308.4  22.7  1.2  305.0  14.5  8.80  100206563  63.9  13.9  4.5  47.2  13.2  16.24  100414097  185.7  15.7  1.0  192.1  8.7  9.38  110106893  155.3  38.6  9.3  134.2  47.0  17.53  110112934  25.9  44.0  4.7  10.6  64.462  22.19  120302080  147.0  25.2  7.7  122.5  29.7  22.15  120624309  0.7  −6.3  1.2  4.8  7.2376  14.14  121128212  278.8  41.6  1.5  300.6  54.3  19.20  GRB  DOL RA  DOL Dec.  DOL Err  True RA  True Dec.  Separation angle  080725435  123.1  −23.1  2.2  121.7  −14.0  9.19  081101491  80.8  11.6  10.1  95.8  −0.1  18.94  090927422  67.6  −67.6  12.1  344.0  −71.0  27.37  100116897  308.4  22.7  1.2  305.0  14.5  8.80  100206563  63.9  13.9  4.5  47.2  13.2  16.24  100414097  185.7  15.7  1.0  192.1  8.7  9.38  110106893  155.3  38.6  9.3  134.2  47.0  17.53  110112934  25.9  44.0  4.7  10.6  64.462  22.19  120302080  147.0  25.2  7.7  122.5  29.7  22.15  120624309  0.7  −6.3  1.2  4.8  7.2376  14.14  121128212  278.8  41.6  1.5  300.6  54.3  19.20  View Large Focusing on the worst case, GRB 121128212, we can specify a few important details about what provides a good location. For the case of using all detectors to fit the spectrum and location, we recover a position similar to the DOL (see Fig. 7). However, when a subset of detectors chosen by starting with the brightest NaI and then selecting detectors via the legal detector set table (see Fig. 6), we find perfect agreement with the true source. Finally, when we mimic the so-called template spectra, we similarly find a position consistent with the DOL. This immediately points out two very important points: obtaining a correct position requires the proper spectral model and not a fixed template; using 12 NaI detectors for spectral and/or location analysis can lead to systematics. Figure 7. View largeDownload slide The location plots from GRB 121128212. From left to right, the results of using all 14 NaI detectors, the result from assuming a template spectrum and the balrog result. The 1σ and 2σ balrog contours are shown in green and purple, respectively. The blue shaded region represents the portion of the sky occulted by the Earth. The various coloured circles are the 60° FOV of the GBM detectors that view the best-fitting balrog position. Figure 7. View largeDownload slide The location plots from GRB 121128212. From left to right, the results of using all 14 NaI detectors, the result from assuming a template spectrum and the balrog result. The 1σ and 2σ balrog contours are shown in green and purple, respectively. The blue shaded region represents the portion of the sky occulted by the Earth. The various coloured circles are the 60° FOV of the GBM detectors that view the best-fitting balrog position. The last point indicates that the calculated off-axis response of the GBM NaI detectors may contain inaccuracies, otherwise, using all 12 NaI detectors would result in an accurate location as shown in our simulations. It is a long-term practice (Goldstein et al. 2012) that for spectral analysis, only a subset of detectors with viewing angles less than 60° should be used due to uncertainties in the GBM NaI off-axis response. Similar attention to detector selection must be taken when fitting for location. Unfortunately, without knowing a location a priori, one must iterate the fitting process until a good match between the detectors viewing the location and those used is found. To demonstrate this, we look at the case of GRB 120624309 whose location that is not within the 60° field of view (FOV) of any NaI, but has detectable signal in all NaI detectors. We iterate through the detector selection and use the value of the marginal likelihood (Z or Bayesian evidence) to choose the best location. We find that certain detector combinations greatly alter the location while others simply change the contours slightly. This can be troublesome in real-time follow-up, and demands better modelling of GBM DRMs to be addressed. In statistical terms, different interval, background and detector selections are different statistical models and highly non-nested. Hence likelihood ratio tests (such as those based on comparing the minimum χ2 of two models) are not valid and should not be used to decide between possible source locations. Still, it should be noted that the use of improper priors such as the uniform priors we assume for spectral parameters makes the marginal likelihood insensitive for model comparison. The performance of the balrog is demonstrated for the entire set of 10 GRBs with the largest offset between the estimated and true positions (see Table 2). We show that we can accurately recover the true positions of all the GRBs within the 95 per cent highest posterior density (HPD) Bayesian credible region. The balrog practically eliminates the systematics in GBM GRB locations. In some cases, the locations were made more accurate via larger error regions, while in many cases, the error region was simply shifted to the true location. In practice, this depends on the intensity of the source. For TRIGDAT data, we find that a CPL model is sufficient for determining a proper location. The medium and soft templates always resulted in a systematic offset unless the true spectral model was similar to the template. In a single case (GRB 120624309), a power-law SED provided the best location that was also indicated by having the highest marginal likelihood. The remaining fits are in shown in Fig. A2 of Appendix A. Table 2. balrog systematic location GRB sample. GRB  True RA  True Dec.  balrog RA  balrog Dec.  balrog separation angle  DOL separation angle  080725435  121.7  −14.0  114.95 : 125.30  −26.62 : −5.93  5.49  9.2  081101491  95.8  −0.1  64.67 : 107.84  −2.63 : 20.78  8.97  18.94  090927422  344.0  −71.0  −82.49 : 110.39  −87.40 : −55.35  22.9  27.38  100116897  305.0  14.5  303.35 : 307.75  13.28 : 16.82  0.6  8.81  100206563  47.2  13.2  34.00 : 51.68  7.08 : 17.12  1.74  16.25  100414097  192.1  8.7  188.61 : 193.40  8.77 : 13.95  3.14  9.38  110112934  10.6  64.462  10.82 : 62.79  51.54 : 72.23  13.74  22.19  120302080  122.5  29.7  4.45 : 218.29  −76.84 : 48.43  30.71  22.16  120624309  4.8  7.2376  −0.43 : 5.61  0.63 : 7.26  3.81  14.14  121128212  300.6  54.3  299.38 : 304.32  51.57 : 54.59  1.38  19.21  GRB  True RA  True Dec.  balrog RA  balrog Dec.  balrog separation angle  DOL separation angle  080725435  121.7  −14.0  114.95 : 125.30  −26.62 : −5.93  5.49  9.2  081101491  95.8  −0.1  64.67 : 107.84  −2.63 : 20.78  8.97  18.94  090927422  344.0  −71.0  −82.49 : 110.39  −87.40 : −55.35  22.9  27.38  100116897  305.0  14.5  303.35 : 307.75  13.28 : 16.82  0.6  8.81  100206563  47.2  13.2  34.00 : 51.68  7.08 : 17.12  1.74  16.25  100414097  192.1  8.7  188.61 : 193.40  8.77 : 13.95  3.14  9.38  110112934  10.6  64.462  10.82 : 62.79  51.54 : 72.23  13.74  22.19  120302080  122.5  29.7  4.45 : 218.29  −76.84 : 48.43  30.71  22.16  120624309  4.8  7.2376  −0.43 : 5.61  0.63 : 7.26  3.81  14.14  121128212  300.6  54.3  299.38 : 304.32  51.57 : 54.59  1.38  19.21  Note. Separation angles for balrog are determined from the maximum likelihood point. View Large 4.2 GRBs well localized by GBM DOL We examine 10 GRBs that are well localized by the DOL, i.e. GRBs with angular separation between the DOL location and the true location smaller than the DOL 1σ error circle radius from Connaughton et al. (2015). Table 3 details the results. We discarded GRB 110625881 because it has no detector with 5σ signal within 60° (see Section 7 for a discussion on our final procedure). The same models used in Section 4.1 are applied to these GRBs as well. Table 3. balrog good location GRB sample. GRB  True RA  True Dec.  balrog RA  balrog Dec.  balrog separation  GBM separation  080723557  176.8  −60.2  162.99 : 180.15  −62.54 : −55.97  2.7  0.98  091003191  251.5  36.6  250.89 : 251.80  36.02 : 36.80  0.25  0.74  100704149  133.6  −24.2  127.04 : 143.25  −25.85 : −19.71  3.11  0.7  100906576  28.7  55.63  7.46 : 39.41  39.42 : 60.38  5.87  0.59  111103441  327.1  −10.5  326.19 : 333.22  −11.57 : 8.30  9.81  0.51  120119170  120.0  −9.1  118.86 : 120.97  −11.02 : −8.56  0.77  0.99  120512112  325.6  13.6  322.10 : 326.62  11.91 : 18.41  1.85  1.4  120709883  318.4  −50.0  305.01 : 335.24  −60.05 : −48.08  4.33  0.96  121125356  228.5  55.3  210.40 : 236.90  44.46 : 58.36  4.49  0.46  130518580  355.7  47.5  353.63 : 356.70  46.10 : 47.68  0.72  0.64  GRB  True RA  True Dec.  balrog RA  balrog Dec.  balrog separation  GBM separation  080723557  176.8  −60.2  162.99 : 180.15  −62.54 : −55.97  2.7  0.98  091003191  251.5  36.6  250.89 : 251.80  36.02 : 36.80  0.25  0.74  100704149  133.6  −24.2  127.04 : 143.25  −25.85 : −19.71  3.11  0.7  100906576  28.7  55.63  7.46 : 39.41  39.42 : 60.38  5.87  0.59  111103441  327.1  −10.5  326.19 : 333.22  −11.57 : 8.30  9.81  0.51  120119170  120.0  −9.1  118.86 : 120.97  −11.02 : −8.56  0.77  0.99  120512112  325.6  13.6  322.10 : 326.62  11.91 : 18.41  1.85  1.4  120709883  318.4  −50.0  305.01 : 335.24  −60.05 : −48.08  4.33  0.96  121125356  228.5  55.3  210.40 : 236.90  44.46 : 58.36  4.49  0.46  130518580  355.7  47.5  353.63 : 356.70  46.10 : 47.68  0.72  0.64  Note. Separation angles for balrog are determined from the maximum likelihood point. View Large Table 3 shows that balrog is able to recover all the 10 well-localized GRBs. We found that for GRB 091003191, balrog is able to recover the true location with subdegree precision. Such accurate (i.e. the error contour covers the true location) and precise (i.e. small error region) gamma-ray localization of GRBs can be used for X-ray and optical afterglow follow-up observations. It is observed that using all 12 NaI detectors in the localization results in less accurate and precise locations. The fact that the CPL gives more accurate locations in the 10 bursts examined in Section. 4.1, as well as with this sample, indicates again that the fixed templates may result in inaccurate locations unless they are close to the intrinsic GRB spectrum. 4.3 Summary of balrog localization Using the model checking framework developed in Section 3.1, we perform model checking of the balrog on real data. Once again, we are in the fortunate situation of knowing the true locations via follow-up X-ray and optical observations of the GRBs in our sample (Connaughton et al. 2015). Fig. 8 compares ρtrue to the distribution of ρsim for both samples. balrog performs equally well on both samples. Still, the model does not perfectly predict the true locations. Fig. 9 shows the angular separation from the true position to the maximum likelihood estimate from balrog and GBM for the systematically off sample and the entire sample. It is important to point out that for the systematically off sample, while the balrog angular separation can be very large, the posterior distribution generally covers the true source position (which is often outside the DOL error radius). Figure 8. View largeDownload slide The model checking for the systematically off (left) and well-located sample (right) from the balrog posterior. Figure 8. View largeDownload slide The model checking for the systematically off (left) and well-located sample (right) from the balrog posterior. Figure 9. View largeDownload slide Left: histogram of balrog and DOL angular separation between the best-fitting location to the true location. Note that all the DOL locations systematically miss the true location, though balrog captures all within the 95 per cent posterior. Right: the same as the left-hand panel, but for all locations. The DOL's tail is longer and includes GRBs that were not located within its error contours. Figure 9. View largeDownload slide Left: histogram of balrog and DOL angular separation between the best-fitting location to the true location. Note that all the DOL locations systematically miss the true location, though balrog captures all within the 95 per cent posterior. Right: the same as the left-hand panel, but for all locations. The DOL's tail is longer and includes GRBs that were not located within its error contours. To further emphasize how balrog compares to the DOL, we perform the same check on the systematically off and good sample using the DOL locations. A major caveat is that the DOL locations accessible via Connaughton et al. (2015) are not the result of Bayesian sampling, and do not possess the non-Gaussian error regions calculated in the DOL. Therefore, we assume the reported location is the centre of a multivariate Gaussian with standard deviation equal to the reported error radius. This forms our pseudo-posterior for each DOL location. We then perform the same model checking analysis as before for the two samples. Fig. 10 shows that the systematically off sample analysed with the DOL performs significantly worse than the balrog at predicting the true locations. Moreover, the tightly peaked shape of the DOL posterior contrasts with the broader distribution of the balrog. This translates to the fact that while the balrog does not perfectly reconstruct the location in all cases (for reasons relating to source geometry, detector calibration, source intensity, etc.), the method admits its ignorance with wider error regions (see Fig. A2). On the other hand, the DOL performs very well for those GRBs that it located accurately as expected. These GRBs are part of the so-called core distribution of the DOL location accuracy distribution as noted in Connaughton et al. (2015). Figure 10. View largeDownload slide The model checking for the systematically off (left) and well-located sample (right) from the assumed DOL posterior. Note the different axis scales. Figure 10. View largeDownload slide The model checking for the systematically off (left) and well-located sample (right) from the assumed DOL posterior. Note the different axis scales. 5 EFFECT ON ESTIMATES OF GRB SPECTRA The uncertainty in location can be viewed as uncertainty in the known DRM for each detector for that position used in a spectral analysis. Just as the background should not be fixed in spectral analysis, so should the uncertainty in the DRM be left free. Therefore, we investigate how this added location freedom can affect the spectral analysis of GRBs. For this study, we use GBM TTE data that have 128-channel resolution. The balrog can produce DRMs for an arbitrary number of channels. We focus on the well-known GRB 080916C (Abdo et al. 2009) that has an external and precisely measured location (RA = 119$$_{.}^{\circ}$$8, Dec. = −56$$_{.}^{\circ}$$6). We break the study into two parts: time-integrated and time-resolved spectral analysis. In both, we focus on detector selection and purported claims of extra spectral components (e.g. Guiriec et al. 2015). While the GRB was also seen by the Fermi Large Area Telescope (LAT), we neglect these data as they are outside the scope of GBM locations and spectra. 5.1 Time-integrated analysis The interval T0+0.-71. s. is chosen for time-integrated analysis (T0 is the GBM trigger time). The detectors with favourable viewing angles are NaI detectors 0, 3, 4, 6, 7 and BGO 0. First we fit the Band function with these detector selections. Alternative SED choices have no effect on correcting the location. We find that the optimal detector selection is NaI detectors 0, 3, 4 and BGO 0 that result in a correct location with large errors (see Table 4). With this detector selection, we additionally fit Band+blackbody (BB) and SBPL law SEDs. We find that all models have the same marginal likelihood and similar locations. Therefore, the simpler Band model is preferred. We note that the Band+BB model found in this analysis differs from that found in Guiriec et al. (2015) with the blackbody's temperature highly unconstrained. This is likely due to the blackbody acting to modify the low-energy slope to accommodate the localization parameters. It has been shown that blackbodies in time-integrated spectra may arise from spectral evolution (Burgess 2014). Yet, this analysis shows that the narrow blackbody spectrum may be absorbing response deficiencies in the GBM DRMs when their uncertainty is not accounted for as is done in the standard GRB spectral analysis. Fig. 11 displays the spectra of all three SED fits. Figure 11. View largeDownload slide A νFν comparison of SEDs used for the integrated spectrum of GRB 080916C. The Band and Band+BB spectra are very similar due to the fact that the blackbody component contributes little to the spectrum. This indicates that the blackbody's contribution in Guiriec et al. (2015) may be due to neglecting the uncertainty of the DRMs and location. Figure 11. View largeDownload slide A νFν comparison of SEDs used for the integrated spectrum of GRB 080916C. The Band and Band+BB spectra are very similar due to the fact that the blackbody component contributes little to the spectrum. This indicates that the blackbody's contribution in Guiriec et al. (2015) may be due to neglecting the uncertainty of the DRMs and location. Table 4. GRB 080916C time-integrated balrog parameters (95 per cent HPD credible interval). Model  RA  Dec.  Ep  α  β  kT  Break scale  log  Z  Banda  132.8 : 136.5  −51.3 : −48.7  276.6 : 370.2  −0.91 : −0.82  −2.24 : −1.94  …  …  −2831.4  Band  120.9 : 144.8  −59.5 : −36.0  296.4 : 416.7  −0.95 : −0.86  −2.28 : −1.94  …  …  −1881.4  Band+BB  119.7 : 145.1  −60.2 : −36.0  238.6 : 422.9  −0.96 : −0.74  −2.26 : −1.92  1.0 : 325.0  …  −1881.1  SBPL  114.8 : 144.3  −62.2 : −37.0  136.2 : 194.3  −1.07 : −1.00  −2.07 : −1.82  …  1.78E-03 : 3.51E-01  −1878.8  Model  RA  Dec.  Ep  α  β  kT  Break scale  log  Z  Banda  132.8 : 136.5  −51.3 : −48.7  276.6 : 370.2  −0.91 : −0.82  −2.24 : −1.94  …  …  −2831.4  Band  120.9 : 144.8  −59.5 : −36.0  296.4 : 416.7  −0.95 : −0.86  −2.28 : −1.94  …  …  −1881.4  Band+BB  119.7 : 145.1  −60.2 : −36.0  238.6 : 422.9  −0.96 : −0.74  −2.26 : −1.92  1.0 : 325.0  …  −1881.1  SBPL  114.8 : 144.3  −62.2 : −37.0  136.2 : 194.3  −1.07 : −1.00  −2.07 : −1.82  …  1.78E-03 : 3.51E-01  −1878.8  Note. aFit including NaI 6 and 7. View Large We now fit the same data assuming no uncertainty in location, i.e. with a fixed DRM corresponding to the known location. Table 5 shows the recovered parameters corresponding to the fixed DRM fits. There are two important observations: the fitted parameters change when using a fixed DRM resulting in softer low-energy slopes and higher Ep and the Band+BB fit results in the well-known kT ∼ 40 keV blackbody with preferential marginal likelihood (Δlog Z ≃ 6) for the model over the Band function. This demonstrates that an overly precise location can induce apparent extra spectral components to compensate. Table 5. GRB 080916C time-integrated fixed DRM parameters (95 per cent HPD credible interval). Model  Ep  α  β  kT  Break scale  log  Z  Band  413.2 : 548.4  −1.02 : −0.95  −2.29 : −2.00  …  …  −237.4  Band+BB  955.9 : 1545.8  −1.28 : −1.16  −2.89: −2.06  41.2 : 48.64  …  −231.2  SBPL  205.1: 249.3  −1.09 : −1.11  −2.13 : −1.88  …  5.21E-03 : 3.98E-01  −234.4  Model  Ep  α  β  kT  Break scale  log  Z  Band  413.2 : 548.4  −1.02 : −0.95  −2.29 : −2.00  …  …  −237.4  Band+BB  955.9 : 1545.8  −1.28 : −1.16  −2.89: −2.06  41.2 : 48.64  …  −231.2  SBPL  205.1: 249.3  −1.09 : −1.11  −2.13 : −1.88  …  5.21E-03 : 3.98E-01  −234.4  View Large 5.2 Time-resolved analysis The high intensity of GRB 080916C caused Fermi to repoint during the burst (Abdo et al. 2009). The balrog assumed one spacecraft orientation when fitting spectra and locations, though we have designed it with the ability to fit multiple spectra simultaneously all sharing the same GRB location. We do not test this capability here and leave it for a future study. However, the 71 s duration used in the time-integrated analysis may introduce systematics as the spacecraft is repointing and a single spacecraft orientation is not valid. We can instead look at time-resolved spectra and use the spacecraft orientation corresponding to individual intervals to analyse the spectra. Ideally, all intervals should give similar locations. Using the same detector selection as was used in the time-integrated analysis, we choose the first eight intervals identified in the intriguing fine-time spectroscopy of Guiriec et al. (2015) and fit them with Band and CPL+BB SEDs. We additionally fit with the new GRB model identified in Guiriec et al. (2015) that consists of three spectral components, a CPL, a blackbody and a power law, where the CPL and PL have their spectral indices fixed to −0.7 and −1.5, respectively. The results from each time interval are detailed in Table 6. As with the time-integrated analysis, there was no significant Bayesian evidence for extra components in the spectra and a Band function fit the spectra sufficiently. The first interval provided an incorrect location as NaI 0 was not pointed close enough to the true location until after Fermi repointed. Therefore, we did not use NaI 0 in the first interval. Interestingly, while the marginal likelihood did not support the three-component model of Guiriec et al. (2015), the model resulted in inaccurate locations in several time intervals (see Fig. 12). Because of its fixed low-energy slope, the model's flexibility to fit the data is passed to the freedom in the location parameters, which adjusts to the poor fitting of the spectral model by giving an incorrect location. The marginal distributions of the Band fit and the CPL+BB+PL (Figs 13 and 14) demonstrate that while for the Band function properly constrained distributions are obtained, the more complex CPL+BB+PL model's parameters fill out their priors and hence lead to lower Bayesian evidence. One could argue that the incredibly small time intervals do not allow for more complex models to be statistically significant, however, the inability of the model to locate the GRB also makes it difficult to explain it as the true spectrum. However, recent work in Guiriec et al. (2016) shows that the model does help to explain the optical data of some GRBs, and hence, further investigation is required. This model is not able to explain the data of GRB 080916C. Figure 12. View largeDownload slide Locations produced from fitting the CPL+BB+PL function (left) and the Band model (right) to the interval T0+3.8–4.3 of GRB 080916C. The CPL+BB+PL model provides an inaccurate location. Figure 12. View largeDownload slide Locations produced from fitting the CPL+BB+PL function (left) and the Band model (right) to the interval T0+3.8–4.3 of GRB 080916C. The CPL+BB+PL model provides an inaccurate location. Figure 13. View largeDownload slide The marginal distribution for the Band (blue) and CPL+BB+PL (red) fit for the time interval T0+3.8–4.3 for GRB 080916C. K represents normalizations for the various spectral components. It is evident that for the CPL+BB+PL model, the CPL component is consistent with zero, allowing the BB to dominate the spectral peak. The differences in the location parameters are also clearly demonstrated. Figure 13. View largeDownload slide The marginal distribution for the Band (blue) and CPL+BB+PL (red) fit for the time interval T0+3.8–4.3 for GRB 080916C. K represents normalizations for the various spectral components. It is evident that for the CPL+BB+PL model, the CPL component is consistent with zero, allowing the BB to dominate the spectral peak. The differences in the location parameters are also clearly demonstrated. Figure 14. View largeDownload slide A νFν comparison between the Band and CPL+BB+PL for the interval T0+3.8–4.3. The width of the curves demonstrates the 95 per cent HPD marginal distributions including the error associated with location and, implicitly, the calibration. Figure 14. View largeDownload slide A νFν comparison between the Band and CPL+BB+PL for the interval T0+3.8–4.3. The width of the curves demonstrates the 95 per cent HPD marginal distributions including the error associated with location and, implicitly, the calibration. Table 6. GRB 080916C time-resolved parameters (95 per cent HPD credible interval). Time from T0  Model  Ep (keV)  α  β  kT (keV)  δ  RA  Dec.  log Z  −0.1 : 0.7  Band  154.6 : 618.9  −0.97 : −0.41  −4.86 : −1.89  …  …  82.5 : 157.9  −76.1 : −9.9  −346.0    CPL+BB  169.5 : 643.2  −0.97 : −0.43  …  1.00 : 314.79  …  89.1 : 158.6  −76.2 : −12.5  −346.4    CPL+BB+PL  22.2 : 97694.0  −0.7 (fix)  …  44.94 : 65.33  −1.5 (fix)  78.6 : 160.0  −75.1 : −7.1  −348.8  0.7 : 1.2  Band  169.5 : 425.5  −0.76 : −0.32  −5.00 : −2.18  …  …  81.4 : 154.9  −71.3 : −18.8  −326.2    CPL+BB  183.9 : 419.7  −0.76 : −0.33  …  1.00 : 217.38  …  81.6 : 155.6  −71.4 : −17.5  −326.5    CPL+BB+PL  16.4 : 82608.5  −0.7 (fix)  …  49.67 : 66.65  −1.5 (fix)  75.7 : 158.9  −71.3 : −6.2  −332.6  1.2 : 1.7  Band  625.8 : 1710.3  −0.97 : −0.71  −5.00 : −2.41  …  …  117.7 : 158.6  −64.7 : −14.3  −348.4    CPL+BB  169.9 : 2347.1  −0.99 : −0.40  …  6.77 : 797.42  …  110.0 : 158.2  −68.1 : −16.1  −347.6    CPL+BB+PL  17.4 : 98087.7  −0.7 (fix)  …  47.23 : 78.55  −1.5 (fix)  156.7 : 164.3  −5.7 : 9.2  −364.7  1.7 : 2.3  Band  93.2 : 255.7  −0.64 : −0.02  −4.20 : −1.77  …  …  92.1 : 155.9  −67.7 : −16.4  −399.3    CPL+BB  119.0 : 378.0  −0.71 : −0.19  …  1.01 : 569.82  …  96.1 : 154.7  −67.1 : −17.8  −400.1    CPL+BB+PL  15.0 : 88705.0  −0.7 (fix)  …  40.66 : 52.89  −1.5 (fix)  86.8 : 160.4  −66.6 : −0.8  −407.0  2.3 : 2.8  Band  270.2 : 634.1  −0.89 : −0.58  −5.00 : −2.45  …  …  100.2 : 154.5  −70.4 : −26.7  −336.0    CPL+BB  126.3 : 1513.9  −1.07 : −0.40  …  1.25 : 373.46  …  105.1 : 154.6  −69.5 : −26.5  −335.6    CPL+BB+PL  10.1 : 53433.7  −0.7 (fix)  …  47.46 : 60.70  −1.5 (fix)  74.7 : 153.8  −72.9 : −20.8  −346.8  2.8 : 3.4  Band  238.9 : 549.0  −0.86 : −0.54  −5.00 : −2.12  …  …  110.3 : 154.9  −66.2 : −21.2  −396.5    CPL+BB  267.9 : 566.2  −0.86 : −0.56  …  1.00 : 263.49  …  110.3 : 154.3  −65.1 : −20.4  −396.8    CPL+BB+PL  11.8 : 67813.5  −0.7 (fix)  …  52.77 : 68.14  −1.5 (fix)  91.7 : 158.5  −67.8 : −9.5  −401.8  3.4 : 3.8  Band  78.8 : 419.6  −0.87 : 0.02  −4.68 : −1.71  …  …  103.5 : 163.0  −66.2 : 2.2  −256.2    CPL+BB  86.2 : 767.5  −0.96 : −0.25  …  1.00 : 353.30  …  113.5 : 163.0  −63.0 : 0.4  −256.6    CPL+BB+PL  12.3 : 62089.7  −0.7 (fix)  …  36.52 : 55.27  −1.5 (fix)  151.8 : 163.0  −11.1 : 8.8  −262.5  3.8 : 4.3  Band  248.5 : 871.2  −0.94 : −0.57  −4.60 : −1.81  …  …  88.0 : 154.3  −68.6 : −20.0  −338.7    CPL+BB  227.9 : 2254.0  −1.02 : −0.57  …  1.22 : 808.97  …  93.1 : 153.9  −68.5 : −21.7  −339.1    CPL+BB+PL  10.5 : 57624.8  −0.7 (fix)  …  47.41 : 75.23  −1.5 (fix)  78.2 : 161.7  −68.0 : 3.6  −348.1  Time from T0  Model  Ep (keV)  α  β  kT (keV)  δ  RA  Dec.  log Z  −0.1 : 0.7  Band  154.6 : 618.9  −0.97 : −0.41  −4.86 : −1.89  …  …  82.5 : 157.9  −76.1 : −9.9  −346.0    CPL+BB  169.5 : 643.2  −0.97 : −0.43  …  1.00 : 314.79  …  89.1 : 158.6  −76.2 : −12.5  −346.4    CPL+BB+PL  22.2 : 97694.0  −0.7 (fix)  …  44.94 : 65.33  −1.5 (fix)  78.6 : 160.0  −75.1 : −7.1  −348.8  0.7 : 1.2  Band  169.5 : 425.5  −0.76 : −0.32  −5.00 : −2.18  …  …  81.4 : 154.9  −71.3 : −18.8  −326.2    CPL+BB  183.9 : 419.7  −0.76 : −0.33  …  1.00 : 217.38  …  81.6 : 155.6  −71.4 : −17.5  −326.5    CPL+BB+PL  16.4 : 82608.5  −0.7 (fix)  …  49.67 : 66.65  −1.5 (fix)  75.7 : 158.9  −71.3 : −6.2  −332.6  1.2 : 1.7  Band  625.8 : 1710.3  −0.97 : −0.71  −5.00 : −2.41  …  …  117.7 : 158.6  −64.7 : −14.3  −348.4    CPL+BB  169.9 : 2347.1  −0.99 : −0.40  …  6.77 : 797.42  …  110.0 : 158.2  −68.1 : −16.1  −347.6    CPL+BB+PL  17.4 : 98087.7  −0.7 (fix)  …  47.23 : 78.55  −1.5 (fix)  156.7 : 164.3  −5.7 : 9.2  −364.7  1.7 : 2.3  Band  93.2 : 255.7  −0.64 : −0.02  −4.20 : −1.77  …  …  92.1 : 155.9  −67.7 : −16.4  −399.3    CPL+BB  119.0 : 378.0  −0.71 : −0.19  …  1.01 : 569.82  …  96.1 : 154.7  −67.1 : −17.8  −400.1    CPL+BB+PL  15.0 : 88705.0  −0.7 (fix)  …  40.66 : 52.89  −1.5 (fix)  86.8 : 160.4  −66.6 : −0.8  −407.0  2.3 : 2.8  Band  270.2 : 634.1  −0.89 : −0.58  −5.00 : −2.45  …  …  100.2 : 154.5  −70.4 : −26.7  −336.0    CPL+BB  126.3 : 1513.9  −1.07 : −0.40  …  1.25 : 373.46  …  105.1 : 154.6  −69.5 : −26.5  −335.6    CPL+BB+PL  10.1 : 53433.7  −0.7 (fix)  …  47.46 : 60.70  −1.5 (fix)  74.7 : 153.8  −72.9 : −20.8  −346.8  2.8 : 3.4  Band  238.9 : 549.0  −0.86 : −0.54  −5.00 : −2.12  …  …  110.3 : 154.9  −66.2 : −21.2  −396.5    CPL+BB  267.9 : 566.2  −0.86 : −0.56  …  1.00 : 263.49  …  110.3 : 154.3  −65.1 : −20.4  −396.8    CPL+BB+PL  11.8 : 67813.5  −0.7 (fix)  …  52.77 : 68.14  −1.5 (fix)  91.7 : 158.5  −67.8 : −9.5  −401.8  3.4 : 3.8  Band  78.8 : 419.6  −0.87 : 0.02  −4.68 : −1.71  …  …  103.5 : 163.0  −66.2 : 2.2  −256.2    CPL+BB  86.2 : 767.5  −0.96 : −0.25  …  1.00 : 353.30  …  113.5 : 163.0  −63.0 : 0.4  −256.6    CPL+BB+PL  12.3 : 62089.7  −0.7 (fix)  …  36.52 : 55.27  −1.5 (fix)  151.8 : 163.0  −11.1 : 8.8  −262.5  3.8 : 4.3  Band  248.5 : 871.2  −0.94 : −0.57  −4.60 : −1.81  …  …  88.0 : 154.3  −68.6 : −20.0  −338.7    CPL+BB  227.9 : 2254.0  −1.02 : −0.57  …  1.22 : 808.97  …  93.1 : 153.9  −68.5 : −21.7  −339.1    CPL+BB+PL  10.5 : 57624.8  −0.7 (fix)  …  47.41 : 75.23  −1.5 (fix)  78.2 : 161.7  −68.0 : 3.6  −348.1  View Large This raises an interesting paradigm for model selection in GRBs, when the true location is known externally, models, such as the new physical models being fit to data (e.g. Burgess et al. 2014; Vurm & Beloborodov 2016), can be rejected if they do not provide proper location parameters. In future work we will address this with the synchrotron model of Burgess et al. (2014). 6 COMMENT ON THE CASE OF GW  150914 GBM The purported detection of electromagnetic radiation by GBM of Laser Interferometer Gravitational-Wave Observatory (LIGO) event GW  150914 (Abbott et al. 2016a,b; Connaughton et al. 2016) was in contrast with upper limits from the INTEGRAL-SPI/Anti-Coincidence Shield (ACS; Savchenko et al. 2016) and an independent analysis of the GBM data (Greiner et al. 2016). Our analysis in Greiner et al. (2016), which concluded that the data were consistent with purely background, differs from that described in Connaughton et al. (2016) in two ways: we used a different statistical model; and, possibly more importantly, we used only two detectors (NaI 5 and BGO 0) rather than all 14. We have found in our analysis herein that the use of all NaI detectors can result in systematically off locations owing to possible deficiencies in the GBM DRMs. Here, we use the balrog with two data sets: NaI 5 and BGO 0; and all 12 NaI and two BGO detectors. We use PL and CPL SED models to match Connaughton et al. (2016). Fig. 15 shows the results of the analysis when using all detectors with both SEDs. Both 1σ and 2σ regions of the GBM location intersect the LIGO arc. Disregarding the fact that using all 14 GBM detectors has the potential to introduce systematics, we can compare the fluences derived from these fits with the upper limits derived in Savchenko et al. (2016). Marginalizing over the location, Fig. 16 shows the marginal fluence distributions over the interval 70–2000 keV for both spectral functions. The PL function does not have a high-energy cut-off and therefore has less believable limits, although the 95 per cent HPD credible region does include the ACS upper limit. The CPL does have a high-energy cut-off and is clearly violating the ACS upper limits. The fact that systematics exist when using all detectors could propagate into the method of Blackburn et al. (2015) that relies on model rates in all GBM detectors built from the DOL templates discussed herein to estimate the significance of transient events in untriggered GBM data. If such systematics exist, and the event has a different spectrum than assumed from the spectral templates, then it is possible that the significance estimated for the GBM event is incorrect. Figure 15. View largeDownload slide The location of the event associated with GW  150914 using all GBM detectors and a PL photon (left) and CPL SED (right). The black points represent positions along the LIGO arc. Using all detectors results in location contours similar to Connaughton et al. (2016), yet it is possible that these are plagued with systematics. Figure 15. View largeDownload slide The location of the event associated with GW  150914 using all GBM detectors and a PL photon (left) and CPL SED (right). The black points represent positions along the LIGO arc. Using all detectors results in location contours similar to Connaughton et al. (2016), yet it is possible that these are plagued with systematics. Figure 16. View largeDownload slide The marginalized fluence distributions resulting from fits using all GBM detectors to the event associated with GW  150914. The red lines indicate the 95 per cent HPD credible interval for each distribution. The distributions are at odds with the ACS upper limits (Savchenko et al. 2016). Figure 16. View largeDownload slide The marginalized fluence distributions resulting from fits using all GBM detectors to the event associated with GW  150914. The red lines indicate the 95 per cent HPD credible interval for each distribution. The distributions are at odds with the ACS upper limits (Savchenko et al. 2016). Fig. 17 shows the results of fits obtained using data from the NaI 5 and BGO 0 detectors. The location is clearly unconstrained. Additionally, the spectrum is unconstrained. These results combined with the non-detection results of Greiner et al. (2016) point to a background (non-astrophysical) origin of the GBM counts data at the time of the GW  150914. In Greiner et al. (2016), we used two different methods to access the spectra with NaI 5 and BGO 0. We used a Poisson–Gaussian likelihood for fixed positions with maximum likelihood estimation and a Bayesian method where the background and source are fit simultaneously. Both methods found that the signal was consistent with background, the method presented here is similar to the maximum likelihood method of Greiner et al. (2016) except that the location parameters are free. Figure 17. View largeDownload slide The location of the event associated with GW  150914 using only NaI 5 and BGO 0 and a PL photon (left) and CPL SED (right). The black points represent positions along the LIGO arc. This subset of detectors has the best chance of representing the true location as they are those with the best viewing angle to the purported source (Connaughton et al. 2016). Figure 17. View largeDownload slide The location of the event associated with GW  150914 using only NaI 5 and BGO 0 and a PL photon (left) and CPL SED (right). The black points represent positions along the LIGO arc. This subset of detectors has the best chance of representing the true location as they are those with the best viewing angle to the purported source (Connaughton et al. 2016). In light of a new study by Veres et al. (2016), where a CPL was fit to the 14 detector GBM data to understand the plausible physical GRB scenarios that could explain the signal, we re-examine the spectral models of the GBM event. The authors use a method, which is not described, to incorporate multiple locations along the LIGO arc and the Monte Carlo maximum likelihood profile of the parameters. In that study, two sets of these Monte Carlos were created. Each set holds either the νFν peak constant as 1 MeV or the low-energy index at α = −0.42. We repeat this analysis to examine the marginal distributions under the balrog scheme. The locations of each fit are consistent with the LIGO arc at the 2σ level (see Fig. 18). We arrive at slightly different marginal distributions (Fig. 19) most likely due to the fact that we fully marginalize over all locations, though it is difficult to compare the methods without a detailed description from Veres et al. (2016). However, both fits violate the ACS upper limits shown in Fig. 20. Figure 18. View largeDownload slide The location of the event associated with GW  150914 using all GBM detectors repeating the analysis of Veres et al. (2016) where we hold Ep = 1 MeV fixed (left) and α = −0.42 fixed (right). The locations are marginally consistent with the LIGO arc at the 2σ level. Figure 18. View largeDownload slide The location of the event associated with GW  150914 using all GBM detectors repeating the analysis of Veres et al. (2016) where we hold Ep = 1 MeV fixed (left) and α = −0.42 fixed (right). The locations are marginally consistent with the LIGO arc at the 2σ level. Figure 19. View largeDownload slide The marginal distributions of the fixed low-energy index (purple) and fixed Ep (green) from the repeated analysis of Veres et al. (2016). The dotted lines represent the fixed values of each fit that show that the two different fits are inconsistent with each other. The fluence (erg cm−2) is integrated over 10–1000 keV. Figure 19. View largeDownload slide The marginal distributions of the fixed low-energy index (purple) and fixed Ep (green) from the repeated analysis of Veres et al. (2016). The dotted lines represent the fixed values of each fit that show that the two different fits are inconsistent with each other. The fluence (erg cm−2) is integrated over 10–1000 keV. Figure 20. View largeDownload slide The fluences from the fits of the CPL models following the analysis of Veres et al. (2016) compared with the ACS upper limit derived in Savchenko et al. (2016). Figure 20. View largeDownload slide The fluences from the fits of the CPL models following the analysis of Veres et al. (2016) compared with the ACS upper limit derived in Savchenko et al. (2016). 7 DISCUSSION Localizing GRBs is the most important function of GBM in the era of multimessenger astronomy. With the possibility of joint LIGO–GBM observations in the near future, GBM localizations have the potential to drastically shrink the broad localization arcs provided by LIGO leading to higher probability that imaging instruments can zero in on the location of a LIGO trigger. With this goal in mind, we have proposed a modern method to localize GBM GRBs that simultaneously fits the spectrum and location of GRBs via Bayesian nested sampling. The technique eliminates the systematics inherent in the standard GBM localization technique reported in Connaughton et al. (2015), though some issues still remain that may be the result of poor spacecraft modelling or improper spectral model choice. Proper spacecraft modelling can be addressed with better simulations but model choice is a problem that plagues GRB spectroscopy for even the very wise cannot see all ends. In attempting to find the true location of GRBs with systematic offsets via balrog, we found that special care must be taken with both model and detector selection. The following methodology was developed. A subset of detectors with 5σ significance above background (Li & Ma 1983) and less than 60° angular separation from the GBM flight software position are initially used and a location is determined using the CPL model. Once a candidate location has been determined, it is checked that the detectors selected view the final location. If not, the selection is modified and repeated. This process is repeated while also iterating model choice until the highest marginal likelihood is found that we then take to correspond to the true location.Proper location depends heavily on detector selection is indicative of deficiencies in the detector response. This could be explained if: (i) the >60° off-axis response of the individual detectors has some calibration uncertainties; (ii) the Bayesian fitting scheme does not perfectly cover the noise-only detectors; (iii) the simulations include deficiencies because they are based on imperfect mass models of the detectors and the spacecraft and the corresponding spacecraft-related scattering. The effect of these deficiencies has been noted already at the pre-launch source survey, when the full Fermi spacecraft with all GBM detectors mounted was irradiated with a radioactive source. The simulations were never able to reproduce the measurements, in particular for detectors that are at large angles relative to the source beam. Since these mismatches did not show up for the calibration measurements of individual detectors (Bissaldi et al. 2009), we presently suspect the above effect (iii) to be dominant. A deeper investigation into these deficiencies may alleviate the need for careful detector selection that currently makes using the balrog difficult in real-time follow-up. Note that this study only analyses a small fraction of the GRBs listed in Connaughton et al. (2015). However, these are the tail cases (those corresponding to the best and worst localizations made by the DOL) and balrog performs well in both. It is conceivable that for the median of the distribution, balrog will perform equally well. We advocate for a deeper study to access this, but feel that the introduction of the concept herein provides the required impetus to modernize GRB location algorithms. It is also worth restating that this method is a superset of the original DOL method. balrog is computationally expensive and fast location estimates will require its implementation on high-performance computing clusters. A standard location on TRIGDAT data takes ∼1–5 min on a 48-core node without hyperthreading. The main computational bottleneck is the DRM generation. In the future, it may be possible to pre-compute a finer base DRM grid, but the multiplicity of spacecraft orientations that require recomputing of the atmospheric scattering will make such a task difficult. A secondary finding of this study is the effect of including the location errors in spectra. The spectral shape of GRB 080916C for both time-integrated and time-resolved analysis favour a simple Band function rather than multiple components. The new model proposed in Guiriec et al. (2015) is not favoured and results in incorrect locations due to its fixed low-energy slope. To reject a spectral model in this way, it is required that a GRB is localized externally by an instrument such as Swift. In the future, optical/radio counterparts will be much easier to identify given the ∼few deg2 error regions as compared to the ∼100 deg2 error regions (where afterglow identification was possible already in about 50 per cent of the cases; Singer et al. 2013, 2015; Lipunov et al. 2016). While these findings focus on a single GRB, they provide impetus to study these effects on a broader sample that will be the subject of a future work. We reiterate that this opens a new path for GRB model selection. While the Band function's flexibility is likely sufficient to allow for proper location, the push for using physical models in GRB spectroscopy that have less flexibility will allow for quick rejection of these models if they result in inaccurate locations. The systematics identified when using all NaI detectors support the claims of Greiner et al. (2016) regarding the purported signal associated with GW  150914 (Connaughton et al. 2016). While including all NaI detectors results in a location that touches the LIGO arc, it also includes most of the Earth. When using a detector selection that corresponds to the position of the LIGO arc, the location contours encompass the entire sky. This, in combination with the lack of signal found in Greiner et al. (2016), places serious constraints on the astrophysical origin of the event. Furthermore, the balrog can be used for future LIGO events as it provides accurate locations for GBM transients that can be employed in optical follow-up strategies. With these findings, we advocate for adopting the balrog for the location of GBM GRBs as it provides realistic and accurate locations suitable for follow-up strategies. Additionally, the shifts in spectral properties found in GRB 080916C warrant an adoption of this scheme for spectral analysis of GBM GRBs but could also be used for similar instruments. We will investigate the GBM GRB catalogues in a forthcoming publication. Such a scheme could be implemented in modern analysis tools such as the Multi-Mission Maximum Likelihood framework (3ML; Vianello et al. 2015) that would allow the balrog to be used to fit data simultaneously with other instruments while still incorporating the uncertainty in location. ACKNOWLEDGEMENTS This work made use of several open source python packages and we are thankful to the authors of matplotlib (Hunter 2007), astropy (Astropy Collaboration 2013), seaborn (Waskom et al. 2016) and pymultinest (Buchner et al. 2014). The authors are thankful to the authors of the original GBM DRM generation code that helped in designing balrog and to Gudlaugur Jóhannesson for insight into making the code more efficient. Additional statistical discussion with Brandon Anderson and Giacomo Vianello helped steer this work in the right direction. We thank Andreas von Kienlin (MPE) and Marc Kippen (Los Alamos) for discussion on the GBM response calibration and simulation. Footnotes 1 GBM and BATSE DRMs differ from classical instrument DRMs because the effective area and response are combined into one matrix. 2 Our balrog tool is completely unrelated, except in name, to the BaLROG project to examine the dynamical influence of bars in galaxies described by Seidel et al. (2015). 3 The actual expressions given in the pgstat section of Arnaud et al. (2015) are −2 times the logarithm of the likelihoods given here. There are also some terms with no S dependence that have been omitted here and the ambiguous ± sign in equation (10) has been made definite. REFERENCES Abbott B. P. et al.  , 2016a, Phys. Rev. Lett. , 116, 061102 CrossRef Search ADS   Abbott B. P. et al.  , 2016b, ApJ , 826, L13 CrossRef Search ADS   Abdo A. A. et al.  , 2009, Science , 323, 1688 CrossRef Search ADS PubMed  Arnaud K., Dorman B., Gordon C., 2015, Xspec User's Guide, v12.9.0, available at: https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html Astropy Collaboration, 2013, A&A , 558, A33 CrossRef Search ADS   Band D., Matteson J., Ford L., 1993, ApJ , 413, 281 CrossRef Search ADS   Bissaldi E. et al.  , 2009, Exp. Astron. , 24, 47 CrossRef Search ADS   Blackburn L., Briggs M. S., Camp J., Christensen N., Connaughton V., Jenke P., Remillard R. A., Veitch J., 2015, ApJS , 217, 8 CrossRef Search ADS   Briggs M. S. et al.  , 1996, ApJ , 459, 40 CrossRef Search ADS   Buchner J. et al.  , 2014, A&A , 564, A125 CrossRef Search ADS   Burgess J. M., 2014, MNRAS , 445, 2589 CrossRef Search ADS   Burgess J. M. et al.  , 2014, ApJ , 784, 17 CrossRef Search ADS   Connaughton V. et al.  , 2015, ApJS , 216, 32 CrossRef Search ADS   Connaughton V. et al.  , 2016, ApJ , 826, L6 CrossRef Search ADS   Feroz F., Hobson M. P., Bridges M., 2009, MNRAS , 398, 1601 CrossRef Search ADS   Gelman A., Carlin J. B., Stern H. S., Dunson D. B., Vehtari A., Rubin D. B., 2014, Bayesian Data Analysis, 3rd edn. CRC Press, Boca Raton, FL Goldstein A. et al.  , 2012, ApJS , 199, 19 CrossRef Search ADS   Greiner J., Burgess J. M., Savchenko V., Yu H. F., 2016, ApJ , 827, L38 CrossRef Search ADS   Guiriec S. et al.  , 2015, ApJ , 807, 148 CrossRef Search ADS   Guiriec S. et al.  , 2016, ApJ , 831, L8 CrossRef Search ADS   Hunter J. D., 2007, Comput. Sci. Eng. , 9, 90 CrossRef Search ADS   Kippen R. M. et al.  , 2007, in Ritz S., Michelson P., Meegan C. A., eds, AIP Conf. Proc. Vol. 921, The First Glast Symposium. Am. Inst. Phys., New York, p. 590 Li T. P., Ma Y. Q., 1983, ApJ , 272, 317 CrossRef Search ADS   Lipunov V. M. et al.  , 2016, MNRAS , 455, 712 CrossRef Search ADS   Loredo T. J., Lamb D. Q., 1992, in Paciesas W. S., Fishman G. J., eds, AIP Conf. Proc. Vol. 265, Gamma-Ray Bursts. Am. Inst. Phys., New York, p. 414 Meegan C. et al.  , 2009, ApJ , 702, 791 CrossRef Search ADS   Paciesas W. S., Pendelton G. N., Wilson R. B., Fishman G. J., Meegan C. A., 1988, BAAS, 726 Savchenko V. et al.  , 2016, ApJ , 820, L36 CrossRef Search ADS   Seidel M. K., Falcón-Barroso J., Martínez-Valpuesta I., Díaz-García S., Laurikainen E., Salo H., Knapen J. H., 2015, MNRAS , 451, 936 CrossRef Search ADS   Singer L. P. et al.  , 2013, ApJ , 776, L34 CrossRef Search ADS   Singer L. P. et al.  , 2015, ApJ , 806, 52 CrossRef Search ADS   Skilling J., 2004, in Fischer R., Preuss R., Toussaint U. V., eds, AIP Conf. Proc. Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering. Am. Inst. Phys., New York, p. 395 van Dyk D. A., Connors A., Kashyap V. L., 2001, ApJ , 548, 224 CrossRef Search ADS   Veres P., Preece R. D., Goldstein A., Mészáros P., Burns E., Connaughton V., 2016, ApJ , 827, L34 CrossRef Search ADS   Vianello G. et al.  , 2015, preprint (arXiv:1507.08343) Vurm I., Beloborodov A. M., 2016, ApJ , 831, 175 CrossRef Search ADS   Waskom M. et al.  , 2016, seaborn: v0.7.0 (January 2016) , available at: http://dx.doi.org/10.5281/zenodo.45133 APPENDIX A: LOCATION PLOTS Figure A1. View largeDownload slide The plots of the remaining nine GRBs well located by the DOL and balrog. Figure A1. View largeDownload slide The plots of the remaining nine GRBs well located by the DOL and balrog. Figure A2. View largeDownload slide The plots of the remaining nine GRBs with large DOL systematic error. Figure A2. View largeDownload slide The plots of the remaining nine GRBs with large DOL systematic error. © 2017 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

Awakening the BALROG: BAyesian Location Reconstruction Of GRBs

Loading next page...
 
/lp/ou_press/awakening-the-balrog-bayesian-location-reconstruction-of-grbs-0gy8Yb0mjn
Publisher
The Royal Astronomical Society
Copyright
© 2017 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/stx2853
Publisher site
See Article on Publisher Site

Abstract

Abstract The accurate spatial location of gamma-ray bursts (GRBs) is crucial for both accurately characterizing their spectra and follow-up observations by other instruments. The Fermi Gamma-ray Burst Monitor (GBM) has the largest field of view for detecting GRBs as it views the entire unocculted sky, but as a non-imaging instrument it relies on the relative count rates observed in each of its 14 detectors to localize transients. Improving its ability to accurately locate GRBs and other transients is vital to the paradigm of multimessenger astronomy, including the electromagnetic follow-up of gravitational wave signals. Here we present the BAyesian Location Reconstruction Of GRBs (balrog) method for localizing and characterizing GBM transients. Our approach eliminates the systematics of previous approaches by simultaneously fitting for the location and spectrum of a source. It also correctly incorporates the uncertainties in the location of a transient into the spectral parameters and produces reliable positional uncertainties for both well-localized sources and those for which the GBM data cannot effectively constrain the position. While computationally expensive, balrog can be implemented to enable quick follow-up of all GBM transient signals. Also, we identify possible response problems that require attention and caution when using standard, public GBM detector response matrices. Finally, we examine the effects of including the uncertainty in location on the spectral parameters of GRB 080916C. We find that spectral parameters change and no extra components are required when these effects are included in contrast to when we use a fixed location. This finding has the potential to alter both the GRB spectral catalogues and the reported spectral composition of some well-known GRBs. methods: data analysis, methods: statistical, gamma-ray burst: general 1 INTRODUCTION The localization of gamma-ray bursts (GRBs) on the sky has played an important role in the history of understanding these extreme objects. The Burst and Transient Source Experiment (BATSE) helped to confirm the cosmological origin of GRBs by determining the isotropic sky distribution of over 3000 GRBs (Briggs et al. 1996). The Fermi Gamma-ray Burst Monitor (GBM), which can monitor the entire sky (with the exception of the portion occulted by the Earth), is currently providing locations for hundreds of GRBs and other transient gamma-ray events every year. With the advent of gravitational wave astronomy and the potential for electromagnetic follow-up of associated events, GBM could play a crucial role in the new era of multimessenger astronomy. Using differential rates in each of its 14 non-imaging detectors, GBM is able to provide approximate localization of transient events that can then be used by imaging observatories to search for counterparts at other wavelengths. Unfortunately, the standard GBM localization procedure, the Daughter Of Location (DOL), produces offsets of 8°–13° from the subsequently confirmed true locations for many GRBs localized by GBM (Connaughton et al. 2015). The exact cause of this systematic is unknown and could be attributed to several factors such as imperfections in the simulated spacecraft response and/or model. Here we hypothesize that the systematics in the standard GBM localization algorithm are a by-product of the assumption that all GRB spectral energy distributions (SEDs) are alike, and of inaccuracies in the detector response matrix1 (DRM) modelling. The DOL uses a set of three template spectra are assumed in order to calculate a spatial grid of putative observed counts on the sky. This grid is then compared with the real observed counts via a χ2-minimization on source position to produce a best-fitting location. The template spectra are built from the ubiquitous Band function (Band, Matteson & Ford 1993) with three sets of fixed spectral parameters that allow for only the amplitude to be varied. Typically, one template is used to locate a GRB, while the other two are reserved for solar flares and soft-gamma repeaters. Since GRBs exhibit a variety of spectral shapes and temporal spectral evolution (e.g. Burgess 2014), the assumption of a fixed set of spectral shapes could influence inferred location of a GRB if its true SED differs from the templates. To explore this possibility we have developed a method where the location and spectrum are inferred simultaneously via Bayesian sampling using a new method: the BAyesian Location Reconstruction Of GRBs (balrog).2 First, we detail our methodology (Section 2) then demonstrate the viability of the method via simulations (Section 3) and real data (Section 4). Finally, we investigate the effect on spectra introduced by including the uncertainty of location where we use GRB 080916C as a case study (Section 5). 2 INFERENCE METHODOLOGY We analyse GBM data by using Bayesian inference to obtain the posterior distribution of the angular position and emission spectrum of a GRB. This requires a model of the GBM instrument (Section 2.1), a GRB source (Section 2.2) and the sky background (Section 2.3), which can then be combined to obtain the full likelihood for the data set (Section 2.4). Combining the likelihood with the parameter priors (Section 2.5), the implied constraints are then obtained by generating samples from the posterior distribution (Section 2.6). This allows for easy marginalization over nuisance parameters, most obviously integrating out the parameters describing the GRB's spectrum to focus on localization. 2.1 GBM instrumental response The on-ground calibration of GBM was described in Bissaldi et al. (2009), which was used as input for designing the GBM DRMs. Any method of localizing and characterizing sources using GBM must account for the complicated nature of its detectors, which have a broad angular response and suffer from energy dispersion (i.e. there is a significant probability that only a fraction of a detected photon's energy is registered), all of which should be embodied in the DRMs. These effects are accounted for using forward folding, described in detail by Vianello et al. (2015), which provides a probabilistic connection between the actual properties of incident photons and the data recorded by the GBM detectors. A critical simplifying approximation that helps make this problem tractable is that the orientations of the Nd = 14 GBM detectors are perfectly known and constant for the duration of the burst being analysed. This means that the rotation angles (represented as quaternions) that define each detector's orientation can be left implicit, although the formalism presented below could be extended to include the changing orientation of the detectors with the satellite's motion. During the on-source period of a GRB each detector is exposed to a radiation field that here is described by the number flux of photons at the detector (i.e. the rate per unit area, per unit energy, per steradian), n(E, ψ), which depends on both photon energy, E, and angular position, ψ, but is assumed to be time independent. The rate at which photons are registered by detector d (with d ∈ {1, 2, …, Nd}) is n(E, ψ) Ad(E, ψ), where Ad(E, ψ) is the effective area of the detector. The angular dependence of the effective area can largely be understood in terms of the geometrical area presented to a putative source, and so it decreases approximately with the cosine of the angle between ψ and nominal detector orientation, although the output of a full numerical simulation is used in practice. Each photon that is registered then deposits an energy ε ≲ E, a process which can be fully described by the conditional distribution Pr(ε|E, ψ, d, detected). The most important aspect of this energy dispersion is that there is a significant probability that ε is considerably lower than E, due to photons Compton scattering then escaping before all their energy is deposited in the scintillation crystal. This aspect of the detector response has an angular dependence, inextricably linking the probability that a given photon is detected with the energy recorded from it. In the GBM analysis formalism photon detection and energy registration are combined into a single function,   \begin{equation} R_d(\varepsilon , E, \psi ) = A_d(E, \psi ) \, {\rm Pr}(\varepsilon | E, \psi , d, {\rm detected}), \end{equation} (1)which is stored as a DRM and fully characterizes each detector. The energy of an individual photon is not recorded per se, but rather each photon is recorded as being in one of Nc = 128 channels that have a nominal correspondence with energy. The rate at which photons are recorded in channel c (with c ∈ {1, 2, …, Nc} ) is   \begin{equation} S_{d,c} = \int _{E_{{\rm min}, {\rm c}}}^{E_{{\rm max}, {\rm c}}} {\rm d}\varepsilon \int _{0}^\infty {\rm d}E \int {\rm d}\Omega \, n(E, \psi ) \, R_d(\varepsilon , E, \psi ), \end{equation} (2)where Emin, c and Emax, c are the minimum and maximum energies for this channel and the inner integral extends over all photon arrival directions, ψ. 2.1.1 DRM generation GBM DRMs are built via the General Response Simulation System (gress) geant4 Monte Carlo code (Kippen et al. 2007) and separated into two components consisting of a direct DRM and an atmospheric scattering DRM of the 14 detectors for each individual GRB position. A grid of direct DRMs created for 272 points on the sky for each of the 12 GBM sodium iodide (NaI) detectors and two bismuth germinate (BGO) detectors is used as the base response library for DRM generation. These responses exist in an energy compressed format that must be decompressed and added to the atmospheric and satellite internal scattering response via a series of transformations that depend on GBM's position and orientation with respect to the Earth (Paciesas et al. 1988). The directional dependence of the response adds in an extra complication, as it is too computationally expensive to calculate Rd(ε, E, ψ) for all values of ψ over the sphere, so this operation is undertaken only for the source positions proposed during the posterior sampling process (Section 2.6). Therefore, we designed efficient DRM generation software based on the original GBM tool that allows for vectorized and parallel generation of DRMs. The final direct DRM is added to the atmospheric scattering DRM to account for scattering of photons off the Earth's atmosphere (Meegan et al. 2009). balrog has the ability to correctly reject source positions that are occulted by the Earth even though these positions are accounted for in the base DRM data base. The area of the sky that is occulted is calculated in real time using the current orbital altitude of the Fermi satellite and finding the horizon line from the spacecraft zenith. But it is more computationally efficient to do this before generating the full detector response, so the DRM is immediately set to Rd(ε, E, ψ) = 0 for occulted directions. 2.2 Source model For an astronomically distant point source (as GRBs are assumed to be) all photons arrive at the detector from the direction of the source, ψs. This allows the number flux defined in Section 2.1 to be factorized as  \begin{equation} n(E, \psi ) = f(E, \,\,{\phi }_{\rm s}) \, \delta _{\rm D}(\psi , \psi _{\rm s}), \end{equation} (3)where $$\,\,{\phi }_{\rm s}$$ is the vector of parameters (normalization, spectral index, etc.) describing the source's SED and δD(ψ, ψs) is the delta function distribution on the sphere centred on ψs. The presence of the delta function means the predicted count rate from the source given in equation (2) simplifies to be   \begin{equation} S_{d,c} (\psi _{\rm s}, \,\,{\phi }_{\rm s}) = \int _{0}^\infty {\rm d}E \, f(E, \,\,{\phi }_{\rm s}) \int _{E_{{\rm min}, {\rm c}}}^{E_{{\rm max}, {\rm c}}} {\rm d}\varepsilon \, R_d(\varepsilon , E, \psi _{\rm s}), \end{equation} (4)where the dependence on the source parameters has been made explicit. The order of integration has been swapped to emphasize the possibility of pre-calculating the inner integral independent of the source emission parameters. The SED models that are used here are Band, cut-off power law (CPL) and smoothly broken power law (SBPL). The SED parameters $$\,\,{\phi }_{\rm s}$$ and the source position ψs represent the full set of GRB parameters to be constrained using the GBM data. 2.3 Background estimation GBM data include an isotropic sky background and an anisotropic Earth albedo background, both of which are time- and energy dependent. Accounting for the combined background rate from these two effects is a necessary step to inferring the properties of the GRBs detected by the instrument. The ideal approach to this problem would be to simultaneously model the background and source (cf. Loredo & Lamb 1992; van Dyk, Connors & Kashyap 2001), although the fact that the background counts are generally fairly high means that it is a good approximation to estimate the background rate separately. The background is taken to be a polynomial in time, and the coefficients are fit to the off-source light curve in each energy channel, with a Poisson likelihood used to ensure valid results even for high-energy channels with low counts. The result of this process is a full joint posterior distribution of the polynomial coefficients. This posterior is then propagated to the on-source period, the result of which is a data-driven estimate of the background rate, $$\hat{B}_{d,c}$$, and its uncertainty, σd, c, in channel c of detector d. Only the former would be used if a simple background subtraction were sufficient; but, as has been shown by Greiner et al. (2016), failure to include the background uncertainty can lead to qualitatively incorrect inferences about the source of interest. The implied posterior distribution for the background is taken to be a truncated normal of the form   \begin{eqnarray} {\rm Pr}(B_{d,c} | \hat{B}_{d,c}, \sigma _{d,c}) & = & \frac{\Theta (B_{d,c}) }{\left[1 + {\rm erf}\left(\hat{B}_{d,c} / 2^{1/2} \, \sigma _{d,c}\right) \right] / 2}\\ & &\times \frac{1}{(2 \pi )^{1/2} \, \sigma _{d,c}}\, {}{\rm e}^{-\left(B_{d,c} - \hat{B}_{d,c}\right)^2 / \left(2 \, \sigma ^2_{d,c} \right)},\nonumber \end{eqnarray} (5)where erf( · ) is the standard error function and the Heaviside step function Θ( · ) imposes the absolute prior constraint that the background rate cannot be negative. 2.4 Likelihood The likelihood, $${\rm Pr}(\lbrace N_{d,c}\rbrace | \psi _{\rm s}, \,\,{\phi }_{\rm s})$$, is the probability, conditional on a given source position ψs and emission model $$\,\,{\phi }_{\rm s}$$, of obtaining the measured counts, {Nd, c}, during the period (of duration T) that the GBM was observing the GRB. As such, it encodes all the relevant information that can be used to localize and characterize the source. Assuming that each incident photon is only registered once (i.e. in a particular channel of a particular detector), the likelihood has the form   \begin{equation} {\rm Pr}(\lbrace N_{d,c}\rbrace | \psi _{\rm s}, \,\,{\phi }_{\rm s}) = \prod _{d = 1}^{N_{\rm d}} \prod _{c = 1}^{N_{\rm c}} {\rm Pr}\left[N_{d,c} | S_{d,c} (\psi _{\rm s}, \,\,{\phi }_{\rm s})\right], \end{equation} (6)where $$S_{d,c} (\psi _{\rm s}, \,\,{\phi }_{\rm s})$$ is the predicted rate at which source photons are registered in channel c of detector d is given in equation (4). 2.4.1 Single channel contribution As the contribution from each channel and detector combination has the same mathematical form, it is simplest to consider this detector/channel combination without reference to detector number d or channel index c, which are hence left implicit in what follows. If the background rate, B, were known precisely, then the probability of N counts being registered during the on-source period would be given by a simple Poisson distribution as   \begin{equation} {\rm Pr}(N | S, B) = \frac{[(S + B)\, T]^{N} \, {\rm e}^{- (S + B)\, T}}{N!}, \end{equation} (7)where S is the rate of photons being registered from the source. This depends on ψs and $$\,\,{\phi }_{\rm s}$$, but it can be considered as the one model parameter of interest here. The fact that the background is not known perfectly must be incorporated, which formally should be done by including B as an extra parameter to be subsequently marginalized out. Fortunately, however, almost all the available information about B comes from the polynomial fitting to the off-source intervals described in Section 2.3, and not the single on-source measurement. It is hence possible to simply average the Poisson distribution given in equation (7) over the plausible values of the background rate. Given the background estimate $$\hat{B}$$ and associated uncertainty σ, and adopting the truncated normal defined in equation (5), the contribution to the likelihood from a single channel is modified to be   \begin{eqnarray} {\rm Pr}(N | S, \hat{B}, \sigma ) & = & \int _0^\infty {\rm d}B \, {\rm Pr}(B | \hat{B}, \sigma )\, {\rm Pr}(N | S, B)\nonumber \\ & \propto & \int _0^\infty {\rm d}B \, (S + B)^{N} {}{\rm e}^{- [(B + S)\, T + (B - \hat{B})^2 / (2 \, \sigma ^2)]},\nonumber \\ \end{eqnarray} (8)where only the terms with that have a contribution from the source parameters have been retained in the second expression. While the above integral could be evaluated numerically, the fact that it will be calculated numerous times during the posterior sampling process motivates using a closed form approximation instead. The option taken here is to evaluate the profile likelihood, in which, for a given value of S, the integral is taken to be the value of the integrand evaluated at its peak (i.e. for the value of the background rate, $$\tilde{B}$$, that maximizes the integrand). This is the approach implemented by Arnaud, Dorman & Gordon (2015) to obtain the pgstat statistic that is an option in XSPEC. This is equivalent3 to assuming that   \begin{eqnarray} {\rm Pr}(N | S, \hat{B}, \sigma ) \propto \left\lbrace \begin{array}{ll}{\rm e}^{-[ (S + \hat{B})\, T - \sigma ^2 \, T^2 / 2 ]}, & \!\!\!\! {\rm if}\ N = 0, \\ (S + \tilde{B})^{N} \, {}{\rm e}^{- [(S + \tilde{B}) \, T + (\tilde{B} - \hat{B})^2 / (2 \, \sigma ^2)]}, & \!\! {\rm if}\ N > 0, \end{array}\right. \nonumber\\ \end{eqnarray} (9)where in the N > 0 case   \begin{eqnarray} \tilde{B} &=& \frac{1}{2} \! \left\lbrace \hat{B} -S -\sigma ^2 T + \! \left[ (S + \sigma ^2 T - \hat{B})^2\right.\right.\nonumber\\ &&\left.\left. -\, 4 (\sigma ^2 T S - \sigma ^2 N - \hat{B} S) \right]^{1/2} \right\rbrace \! . \end{eqnarray} (10) Independent of whether numerical integration or the profile likelihood is used to evaluate $${\rm Pr}(N | S, \hat{B}, \sigma )$$, there is an increase in the range of values of S for which the likelihood is appreciable. This can lead to the apparently undesirable result that the parameter constraints are similarly broadened, but in fact this is simply a correct reflection of the range of source positions that are consistent with the data and hence which should be considered as potential GRB locations. 2.5 Parameter priors A necessary component of Bayesian parameter estimation is the specification of priors for all parameters. The joint prior in ψs and $$\,\,{\phi }_{\rm s}$$ can reasonably be assumed to factorize (i.e. the plausible SED parameters are independent of position on the sky), so that $${\rm Pr}(\psi _{\rm s}, \,\,{\phi }_{\rm s} | {\rm GRB}) = {\rm Pr}(\psi _{\rm s} | {\rm GRB}) \, {\rm Pr}(\,\,{\phi }_{\rm s} | {\rm GRB})$$. For ψs we assume a prior that is uniform over the celestial sphere, so that Pr(ψs|GRB) = 1/(4π), although there is no operational need to normalize this distribution correctly. For $$\,\,{\phi }_{\rm s}$$ we assume uniform improper priors with appropriate bounds determined for the spectral model at hand (e.g. to ensure that the SED is not negative for positive E). The choice of prior can of course be more physically motivated, however, for our purposes these choices are sufficient because, as we show in Section 4, the marginalized positional posteriors do not have significant prior dependence. 2.6 Posterior distribution and sampling Given the likelihood described in Section 2.4 and the parameter priors given in Section 2.5, the desired posterior distribution has the form   \begin{equation} {\rm Pr}(\psi _{\rm s}, \,\,{\phi }_{\rm s} | \lbrace N_{d,c}\rbrace , {\rm GRB}) \propto {\rm Pr}(\psi _{\rm s}, \,\,{\phi }_{\rm s} | {\rm GRB}) \, {\rm Pr}(\lbrace N_{d,c}\rbrace | \psi _{\rm s}, \,\,{\phi }_{\rm s}).\!\!\! \end{equation} (11)There is no analytic way of determining the combinations of parameter values for which the posterior is high, nor of obtaining the normalizing constant that is needed to marginalize over $$\,\,{\phi }_{\rm s}$$. As the posterior can be multimodal, we use the form of nested sampling (Skilling 2004) as implemented in multinest (Feroz, Hobson & Bridges 2009) to obtain samples from the posterior. A particular utility of representing the posterior distribution as a set of samples is that marginalization is straightforward. In our case, depending on the situation, either $$\,\,{\phi }_{\rm s}$$ or ψs are nuisance parameters and marginalizing over one or other allows us to examine the marginal distribution of the parameter(s) of interest, Pr(ψs|{Nd, c}, GRB) or $${\rm Pr}(\,\,{\phi }_{\rm s} | \lbrace N_{d,c}\rbrace , {\rm GRB})$$, respectively. In order to be useful for timely follow-up, the balrog must be able to localize a transient swiftly. The current implementation can be run on GBM data types including the quick-look TRIGDAT that consists of light curves for each of the 14 GBM detectors each with the nominal 128-channel data rebinned to 8-channel energy resolution, as well as GBM time-tagged event (TTE) data that consist of photon time tags with 128-channel energy resolution (Meegan et al. 2009). Localization on a multicore desktop workstation takes ∼5–10 min for TRIGDAT data and ∼10–20 min for TTE data due to its finer energy resolution. Faster speeds can be achieved with a computing cluster-based installation that is used in this study. We discuss the differences between results with TRIGDAT and TTE data in the following sections. 3 SYNTHETIC LOCATION ANALYSIS To assess the validity of our method, we generate synthetic GRBs with known location and spectrum on a 25-point grid on the sky for a given spacecraft orientation. The CPL SED is used for the synthetic photon spectrum with power-law index of −1 and νFν peak of 300 keV. For each point in the grid, a DRM is generated for the given sky position and the SED is convolved with the DRM and given Poisson variability to produce an 8-channel synthetic data spectrum resembling TRIGDAT data. A simple Poisson background is added to the data obeying a power-law energy spectrum. Using balrog, we fit each simulated GRB in the grid with the CPL model to determine the location and spectrum using all 12 NaI detectors. Fig. 1 shows the distribution of angular separations between the true simulated position and the best-fitting > position. This is always less than 1°and has a mean of ∼0$$_{.}^{\circ}$$1. For real sources, the error could be much larger and depends on the intensity of the GRB. These results depend on knowing the true model of the spectrum and the calibration of the NaI detectors being perfect. As we will see in the following sections, errors in the response model will lead to systematic deviations from the true location when all 12 NaI detectors are used. Figure 1. View largeDownload slide Distribution of angular separations between the true simulated positions and the best-fitting ? position using the CPL SED for the 25-point simulated grid. Figure 1. View largeDownload slide Distribution of angular separations between the true simulated positions and the best-fitting ? position using the CPL SED for the 25-point simulated grid. In order to demonstrate the effect of using fixed templates to determine locations, we fit the simulated grid with a Band function with spectral parameters Ep = 1 MeV, α = 0 and β = −1.5 (i.e. the hard template). The amplitude and location are the only parameters allowed to vary in the fit. Fig. 2 shows that using the hard template introduces a systematic offset in the location of ∼13°, which is comparable to that identified by Connaughton et al. (2015). Moreover, we can examine the angular separation between the best-fitting positions and the location posterior (essentially the error region of the location) for one point on the grid and compare it with the angular separation between the true position and the location posterior. Fig. 3 shows that, when fitting with the CPL model, these two distributions are indistinguishable, which indicates that the error region of the fit is encompassing the true uncertainty in the data. However, while the error region obtained from the hard template is very precise, it is not accurate (i.e. does not encompass the true error region). Figure 2. View largeDownload slide Distribution of angular separations from the simulated positions to the best-fitting 4 position using the hard template SED for the 25-point simulated grid. Figure 2. View largeDownload slide Distribution of angular separations from the simulated positions to the best-fitting 4 position using the hard template SED for the 25-point simulated grid. Figure 3. View largeDownload slide Comparison between distributions of the angular separation of the location posterior and the best-fitting location (blue) and true location (green) for the CPL and hard template SEDs. Figure 3. View largeDownload slide Comparison between distributions of the angular separation of the location posterior and the best-fitting location (blue) and true location (green) for the CPL and hard template SEDs. We can look at how the differences in the locations found by the CPL and hard template models affect the DRM used for spectral analysis. As a toy example, we propagate the location posterior through the DRM generator to create a probabilistic DRM for both models. Fig. 4 illustrates the difference between these two DRMs. The main differences are at high and low energy, which can cause different slopes and hence different spectral parameters to be recovered. However, in the balrog scheme, the DRM is no longer a fixed quantity, but rather a probabilistic quantity weighted by the uncertainty in location. Figure 4. View largeDownload slide A probabilistic difference between the DRMs generated from the location posterior suing the hard template and the CPL SED. The colour scale indicates the difference in effective area in arbitrary units. This serves to illustrate the differences that can arise in the normal public DRMs when an improper model is used to generate DRMs. Figure 4. View largeDownload slide A probabilistic difference between the DRMs generated from the location posterior suing the hard template and the CPL SED. The colour scale indicates the difference in effective area in arbitrary units. This serves to illustrate the differences that can arise in the normal public DRMs when an improper model is used to generate DRMs. 3.1 Posterior checking To assess whether balrog provides GRBs position estimates of GRBs with realistic uncertainties, we need to perform model checking. The standard Bayesian approach is a posterior predictive check (PPC; see e.g. Gelman et al. 2014), in which simulated data sets for the GRB, drawn from the posterior distribution, are compared to the actual data on that source. Here, however, we are in the unusual situation of knowing the true locations of the GRBs from either because there accurate positions from other instruments such as Swift or because the data were simulated with known parameters. That allows us to bypass the data-simulation step and simply assess whether our posterior distributions are consistent with the GRBs’ actual positions. For a single GRB, with known sky position ψtrue, the key statistical quantity is the posterior density at this location, ρtrue = Pr(ψs = ψtrue|{Nd, c}). In cases where the posterior is singly peaked and symmetric this reduces to a measure based on the angular separation between the peak of the posterior and the true position, but is also valid if the posterior has significant degeneracies or is multimodal. The numerical value of ρtrue is not meaningful in isolation, as it is a probability density that is determined by e.g. the scale of the positional uncertainties. It must instead be evaluated in relation to the distribution of ρ values that would be obtained by drawing random positions from the posterior. A sample from this distribution can be obtained by drawing a position ψsim from the posterior distribution Pr(ψs|{Nd, c}) and then evaluating ρsim = Pr(ψs = ψsim|{Nd, c}) (i.e. the density at this simulated position). Repeating this process yields a set of samples from the distribution of ρ values from which ρtrue ought to have been drawn if the model were correct. If the model is badly wrong, then the true position will lie outside the region identified by the posterior and ρtrue will be lower than (almost all of) the simulated values ρsim. Conversely, if the uncertainties are overestimated, then ρtrue will be higher than the simulated values. This comparison can be made on a source-by-source basis, but for a set of NGRB GRBs the ρ values can be multiplied to obtain   \begin{equation} \rho _{\rm true} \equiv \prod _{i = 1}^{N_{\rm GRB}} \rho _{{\rm true},i}, \end{equation} (12)where ρi ≡ Pr(ψs = ψtrue, i|{Nd, c}, i) is evaluated from the posterior distribution for the ith GRB. Similarly, with ρsim, i, j being the jth simulated ρ value drawn from the posterior of the ith GRB, a sample from the expected distribution is given by the product   \begin{equation} \rho _{{\rm sim},j} \equiv \prod _{i = 1}^{N_{\rm GRB}} \rho _{{\rm sim},i,j}. \end{equation} (13) Applying this framework to the simulations, we obtain Fig. 5. As expected from simulations assuming an ideal detector response, the model is accurate at predicting future predictions (ρtrue is not too far to the left) without overfitting the data with the increased uncertainty of the balrog model (ρtrue is not too far to the right). Figure 5. View largeDownload slide The model checking distribution for the simulated GRB sets using the full approach of the balrog. Figure 5. View largeDownload slide The model checking distribution for the simulated GRB sets using the full approach of the balrog. 4 ANALYSIS OF REAL GRBS WITH KNOWN LOCATIONS Accessing the validity of our method depends strongly on our ability to precisely and accurately locate real GRBs that have been located by imaging instruments. We use the table of known GRB locations provided in Connaughton et al. (2015) that additionally allows us to compare our method to GRBs localized by the standard GBM method. We divide the analysis into two sections: GRBs well localized by GBM DOL and GRBs localized with a large systematic deviation by GBM DOL. First, we note a few details about fitting real data. In the above simulations, we assumed the DRMs are perfect and expect to recover the true location and spectral parameters. If there are calibration/simulation effects not properly modelled by the DRMs, then analysis of real data can produce suboptimal results leading to offsets in location and spectral parameters. Similarly, if the spectral model we choose to fit the GRBs does not correspond to the true spectral model of the source, we can expect the location to be affected. In real-time trigger situations, the user of balrog will not know the location a priori. As we will demonstrate in the following sections, the user must select a subset of detectors to obtain correct locations of sources. The GBM team has assembled a table of so-called legal detector sets to aid in choosing proper subsets of detectors for analysis (see Fig. 6; Meegan, team internal information). It will often be required to perform several location runs, with different spectral models and detector subsets to find the best location (see Section 7). balrog provides the marginal likelihood or evidence that can be used to compare between these iterations and choose the best location. It is currently unlikely that an automated system can perform these actions due to the complicated interaction between background, source interval selection, model choice and detector selection. Figure 6. View largeDownload slide Chart of GBM NaI detectors that see the same source based on a simulation of 10  000 GRBs made by the Fermi GBM team. Colour and number correspond to the number of coincident detections between a pair of detectors. It is intended to aid in choosing appropriate sets of detectors for spectral and location analysis. Figure 6. View largeDownload slide Chart of GBM NaI detectors that see the same source based on a simulation of 10  000 GRBs made by the Fermi GBM team. Colour and number correspond to the number of coincident detections between a pair of detectors. It is intended to aid in choosing appropriate sets of detectors for spectral and location analysis. 4.1 GRBs localized with GBM DOL resulting in a systematic deviation We examine 10 GRBs with systematic deviation in the DOL location from the true location to test if balrog can improve upon the result or provide appropriate location contours. We define a systematic deviation such that the angular separation between the DOL location and the true location exceeds the reported DOL 1σ error. We choose the 10 worst cases from Connaughton et al. (2015) that are listed in Table 1. For each of these GRBs, we try three spectral models (Band, CPL and SBPL) and two detector selections: all detectors and a subset of detectors. Table 1. Systematic location GRB sample. GRB  DOL RA  DOL Dec.  DOL Err  True RA  True Dec.  Separation angle  080725435  123.1  −23.1  2.2  121.7  −14.0  9.19  081101491  80.8  11.6  10.1  95.8  −0.1  18.94  090927422  67.6  −67.6  12.1  344.0  −71.0  27.37  100116897  308.4  22.7  1.2  305.0  14.5  8.80  100206563  63.9  13.9  4.5  47.2  13.2  16.24  100414097  185.7  15.7  1.0  192.1  8.7  9.38  110106893  155.3  38.6  9.3  134.2  47.0  17.53  110112934  25.9  44.0  4.7  10.6  64.462  22.19  120302080  147.0  25.2  7.7  122.5  29.7  22.15  120624309  0.7  −6.3  1.2  4.8  7.2376  14.14  121128212  278.8  41.6  1.5  300.6  54.3  19.20  GRB  DOL RA  DOL Dec.  DOL Err  True RA  True Dec.  Separation angle  080725435  123.1  −23.1  2.2  121.7  −14.0  9.19  081101491  80.8  11.6  10.1  95.8  −0.1  18.94  090927422  67.6  −67.6  12.1  344.0  −71.0  27.37  100116897  308.4  22.7  1.2  305.0  14.5  8.80  100206563  63.9  13.9  4.5  47.2  13.2  16.24  100414097  185.7  15.7  1.0  192.1  8.7  9.38  110106893  155.3  38.6  9.3  134.2  47.0  17.53  110112934  25.9  44.0  4.7  10.6  64.462  22.19  120302080  147.0  25.2  7.7  122.5  29.7  22.15  120624309  0.7  −6.3  1.2  4.8  7.2376  14.14  121128212  278.8  41.6  1.5  300.6  54.3  19.20  View Large Focusing on the worst case, GRB 121128212, we can specify a few important details about what provides a good location. For the case of using all detectors to fit the spectrum and location, we recover a position similar to the DOL (see Fig. 7). However, when a subset of detectors chosen by starting with the brightest NaI and then selecting detectors via the legal detector set table (see Fig. 6), we find perfect agreement with the true source. Finally, when we mimic the so-called template spectra, we similarly find a position consistent with the DOL. This immediately points out two very important points: obtaining a correct position requires the proper spectral model and not a fixed template; using 12 NaI detectors for spectral and/or location analysis can lead to systematics. Figure 7. View largeDownload slide The location plots from GRB 121128212. From left to right, the results of using all 14 NaI detectors, the result from assuming a template spectrum and the balrog result. The 1σ and 2σ balrog contours are shown in green and purple, respectively. The blue shaded region represents the portion of the sky occulted by the Earth. The various coloured circles are the 60° FOV of the GBM detectors that view the best-fitting balrog position. Figure 7. View largeDownload slide The location plots from GRB 121128212. From left to right, the results of using all 14 NaI detectors, the result from assuming a template spectrum and the balrog result. The 1σ and 2σ balrog contours are shown in green and purple, respectively. The blue shaded region represents the portion of the sky occulted by the Earth. The various coloured circles are the 60° FOV of the GBM detectors that view the best-fitting balrog position. The last point indicates that the calculated off-axis response of the GBM NaI detectors may contain inaccuracies, otherwise, using all 12 NaI detectors would result in an accurate location as shown in our simulations. It is a long-term practice (Goldstein et al. 2012) that for spectral analysis, only a subset of detectors with viewing angles less than 60° should be used due to uncertainties in the GBM NaI off-axis response. Similar attention to detector selection must be taken when fitting for location. Unfortunately, without knowing a location a priori, one must iterate the fitting process until a good match between the detectors viewing the location and those used is found. To demonstrate this, we look at the case of GRB 120624309 whose location that is not within the 60° field of view (FOV) of any NaI, but has detectable signal in all NaI detectors. We iterate through the detector selection and use the value of the marginal likelihood (Z or Bayesian evidence) to choose the best location. We find that certain detector combinations greatly alter the location while others simply change the contours slightly. This can be troublesome in real-time follow-up, and demands better modelling of GBM DRMs to be addressed. In statistical terms, different interval, background and detector selections are different statistical models and highly non-nested. Hence likelihood ratio tests (such as those based on comparing the minimum χ2 of two models) are not valid and should not be used to decide between possible source locations. Still, it should be noted that the use of improper priors such as the uniform priors we assume for spectral parameters makes the marginal likelihood insensitive for model comparison. The performance of the balrog is demonstrated for the entire set of 10 GRBs with the largest offset between the estimated and true positions (see Table 2). We show that we can accurately recover the true positions of all the GRBs within the 95 per cent highest posterior density (HPD) Bayesian credible region. The balrog practically eliminates the systematics in GBM GRB locations. In some cases, the locations were made more accurate via larger error regions, while in many cases, the error region was simply shifted to the true location. In practice, this depends on the intensity of the source. For TRIGDAT data, we find that a CPL model is sufficient for determining a proper location. The medium and soft templates always resulted in a systematic offset unless the true spectral model was similar to the template. In a single case (GRB 120624309), a power-law SED provided the best location that was also indicated by having the highest marginal likelihood. The remaining fits are in shown in Fig. A2 of Appendix A. Table 2. balrog systematic location GRB sample. GRB  True RA  True Dec.  balrog RA  balrog Dec.  balrog separation angle  DOL separation angle  080725435  121.7  −14.0  114.95 : 125.30  −26.62 : −5.93  5.49  9.2  081101491  95.8  −0.1  64.67 : 107.84  −2.63 : 20.78  8.97  18.94  090927422  344.0  −71.0  −82.49 : 110.39  −87.40 : −55.35  22.9  27.38  100116897  305.0  14.5  303.35 : 307.75  13.28 : 16.82  0.6  8.81  100206563  47.2  13.2  34.00 : 51.68  7.08 : 17.12  1.74  16.25  100414097  192.1  8.7  188.61 : 193.40  8.77 : 13.95  3.14  9.38  110112934  10.6  64.462  10.82 : 62.79  51.54 : 72.23  13.74  22.19  120302080  122.5  29.7  4.45 : 218.29  −76.84 : 48.43  30.71  22.16  120624309  4.8  7.2376  −0.43 : 5.61  0.63 : 7.26  3.81  14.14  121128212  300.6  54.3  299.38 : 304.32  51.57 : 54.59  1.38  19.21  GRB  True RA  True Dec.  balrog RA  balrog Dec.  balrog separation angle  DOL separation angle  080725435  121.7  −14.0  114.95 : 125.30  −26.62 : −5.93  5.49  9.2  081101491  95.8  −0.1  64.67 : 107.84  −2.63 : 20.78  8.97  18.94  090927422  344.0  −71.0  −82.49 : 110.39  −87.40 : −55.35  22.9  27.38  100116897  305.0  14.5  303.35 : 307.75  13.28 : 16.82  0.6  8.81  100206563  47.2  13.2  34.00 : 51.68  7.08 : 17.12  1.74  16.25  100414097  192.1  8.7  188.61 : 193.40  8.77 : 13.95  3.14  9.38  110112934  10.6  64.462  10.82 : 62.79  51.54 : 72.23  13.74  22.19  120302080  122.5  29.7  4.45 : 218.29  −76.84 : 48.43  30.71  22.16  120624309  4.8  7.2376  −0.43 : 5.61  0.63 : 7.26  3.81  14.14  121128212  300.6  54.3  299.38 : 304.32  51.57 : 54.59  1.38  19.21  Note. Separation angles for balrog are determined from the maximum likelihood point. View Large 4.2 GRBs well localized by GBM DOL We examine 10 GRBs that are well localized by the DOL, i.e. GRBs with angular separation between the DOL location and the true location smaller than the DOL 1σ error circle radius from Connaughton et al. (2015). Table 3 details the results. We discarded GRB 110625881 because it has no detector with 5σ signal within 60° (see Section 7 for a discussion on our final procedure). The same models used in Section 4.1 are applied to these GRBs as well. Table 3. balrog good location GRB sample. GRB  True RA  True Dec.  balrog RA  balrog Dec.  balrog separation  GBM separation  080723557  176.8  −60.2  162.99 : 180.15  −62.54 : −55.97  2.7  0.98  091003191  251.5  36.6  250.89 : 251.80  36.02 : 36.80  0.25  0.74  100704149  133.6  −24.2  127.04 : 143.25  −25.85 : −19.71  3.11  0.7  100906576  28.7  55.63  7.46 : 39.41  39.42 : 60.38  5.87  0.59  111103441  327.1  −10.5  326.19 : 333.22  −11.57 : 8.30  9.81  0.51  120119170  120.0  −9.1  118.86 : 120.97  −11.02 : −8.56  0.77  0.99  120512112  325.6  13.6  322.10 : 326.62  11.91 : 18.41  1.85  1.4  120709883  318.4  −50.0  305.01 : 335.24  −60.05 : −48.08  4.33  0.96  121125356  228.5  55.3  210.40 : 236.90  44.46 : 58.36  4.49  0.46  130518580  355.7  47.5  353.63 : 356.70  46.10 : 47.68  0.72  0.64  GRB  True RA  True Dec.  balrog RA  balrog Dec.  balrog separation  GBM separation  080723557  176.8  −60.2  162.99 : 180.15  −62.54 : −55.97  2.7  0.98  091003191  251.5  36.6  250.89 : 251.80  36.02 : 36.80  0.25  0.74  100704149  133.6  −24.2  127.04 : 143.25  −25.85 : −19.71  3.11  0.7  100906576  28.7  55.63  7.46 : 39.41  39.42 : 60.38  5.87  0.59  111103441  327.1  −10.5  326.19 : 333.22  −11.57 : 8.30  9.81  0.51  120119170  120.0  −9.1  118.86 : 120.97  −11.02 : −8.56  0.77  0.99  120512112  325.6  13.6  322.10 : 326.62  11.91 : 18.41  1.85  1.4  120709883  318.4  −50.0  305.01 : 335.24  −60.05 : −48.08  4.33  0.96  121125356  228.5  55.3  210.40 : 236.90  44.46 : 58.36  4.49  0.46  130518580  355.7  47.5  353.63 : 356.70  46.10 : 47.68  0.72  0.64  Note. Separation angles for balrog are determined from the maximum likelihood point. View Large Table 3 shows that balrog is able to recover all the 10 well-localized GRBs. We found that for GRB 091003191, balrog is able to recover the true location with subdegree precision. Such accurate (i.e. the error contour covers the true location) and precise (i.e. small error region) gamma-ray localization of GRBs can be used for X-ray and optical afterglow follow-up observations. It is observed that using all 12 NaI detectors in the localization results in less accurate and precise locations. The fact that the CPL gives more accurate locations in the 10 bursts examined in Section. 4.1, as well as with this sample, indicates again that the fixed templates may result in inaccurate locations unless they are close to the intrinsic GRB spectrum. 4.3 Summary of balrog localization Using the model checking framework developed in Section 3.1, we perform model checking of the balrog on real data. Once again, we are in the fortunate situation of knowing the true locations via follow-up X-ray and optical observations of the GRBs in our sample (Connaughton et al. 2015). Fig. 8 compares ρtrue to the distribution of ρsim for both samples. balrog performs equally well on both samples. Still, the model does not perfectly predict the true locations. Fig. 9 shows the angular separation from the true position to the maximum likelihood estimate from balrog and GBM for the systematically off sample and the entire sample. It is important to point out that for the systematically off sample, while the balrog angular separation can be very large, the posterior distribution generally covers the true source position (which is often outside the DOL error radius). Figure 8. View largeDownload slide The model checking for the systematically off (left) and well-located sample (right) from the balrog posterior. Figure 8. View largeDownload slide The model checking for the systematically off (left) and well-located sample (right) from the balrog posterior. Figure 9. View largeDownload slide Left: histogram of balrog and DOL angular separation between the best-fitting location to the true location. Note that all the DOL locations systematically miss the true location, though balrog captures all within the 95 per cent posterior. Right: the same as the left-hand panel, but for all locations. The DOL's tail is longer and includes GRBs that were not located within its error contours. Figure 9. View largeDownload slide Left: histogram of balrog and DOL angular separation between the best-fitting location to the true location. Note that all the DOL locations systematically miss the true location, though balrog captures all within the 95 per cent posterior. Right: the same as the left-hand panel, but for all locations. The DOL's tail is longer and includes GRBs that were not located within its error contours. To further emphasize how balrog compares to the DOL, we perform the same check on the systematically off and good sample using the DOL locations. A major caveat is that the DOL locations accessible via Connaughton et al. (2015) are not the result of Bayesian sampling, and do not possess the non-Gaussian error regions calculated in the DOL. Therefore, we assume the reported location is the centre of a multivariate Gaussian with standard deviation equal to the reported error radius. This forms our pseudo-posterior for each DOL location. We then perform the same model checking analysis as before for the two samples. Fig. 10 shows that the systematically off sample analysed with the DOL performs significantly worse than the balrog at predicting the true locations. Moreover, the tightly peaked shape of the DOL posterior contrasts with the broader distribution of the balrog. This translates to the fact that while the balrog does not perfectly reconstruct the location in all cases (for reasons relating to source geometry, detector calibration, source intensity, etc.), the method admits its ignorance with wider error regions (see Fig. A2). On the other hand, the DOL performs very well for those GRBs that it located accurately as expected. These GRBs are part of the so-called core distribution of the DOL location accuracy distribution as noted in Connaughton et al. (2015). Figure 10. View largeDownload slide The model checking for the systematically off (left) and well-located sample (right) from the assumed DOL posterior. Note the different axis scales. Figure 10. View largeDownload slide The model checking for the systematically off (left) and well-located sample (right) from the assumed DOL posterior. Note the different axis scales. 5 EFFECT ON ESTIMATES OF GRB SPECTRA The uncertainty in location can be viewed as uncertainty in the known DRM for each detector for that position used in a spectral analysis. Just as the background should not be fixed in spectral analysis, so should the uncertainty in the DRM be left free. Therefore, we investigate how this added location freedom can affect the spectral analysis of GRBs. For this study, we use GBM TTE data that have 128-channel resolution. The balrog can produce DRMs for an arbitrary number of channels. We focus on the well-known GRB 080916C (Abdo et al. 2009) that has an external and precisely measured location (RA = 119$$_{.}^{\circ}$$8, Dec. = −56$$_{.}^{\circ}$$6). We break the study into two parts: time-integrated and time-resolved spectral analysis. In both, we focus on detector selection and purported claims of extra spectral components (e.g. Guiriec et al. 2015). While the GRB was also seen by the Fermi Large Area Telescope (LAT), we neglect these data as they are outside the scope of GBM locations and spectra. 5.1 Time-integrated analysis The interval T0+0.-71. s. is chosen for time-integrated analysis (T0 is the GBM trigger time). The detectors with favourable viewing angles are NaI detectors 0, 3, 4, 6, 7 and BGO 0. First we fit the Band function with these detector selections. Alternative SED choices have no effect on correcting the location. We find that the optimal detector selection is NaI detectors 0, 3, 4 and BGO 0 that result in a correct location with large errors (see Table 4). With this detector selection, we additionally fit Band+blackbody (BB) and SBPL law SEDs. We find that all models have the same marginal likelihood and similar locations. Therefore, the simpler Band model is preferred. We note that the Band+BB model found in this analysis differs from that found in Guiriec et al. (2015) with the blackbody's temperature highly unconstrained. This is likely due to the blackbody acting to modify the low-energy slope to accommodate the localization parameters. It has been shown that blackbodies in time-integrated spectra may arise from spectral evolution (Burgess 2014). Yet, this analysis shows that the narrow blackbody spectrum may be absorbing response deficiencies in the GBM DRMs when their uncertainty is not accounted for as is done in the standard GRB spectral analysis. Fig. 11 displays the spectra of all three SED fits. Figure 11. View largeDownload slide A νFν comparison of SEDs used for the integrated spectrum of GRB 080916C. The Band and Band+BB spectra are very similar due to the fact that the blackbody component contributes little to the spectrum. This indicates that the blackbody's contribution in Guiriec et al. (2015) may be due to neglecting the uncertainty of the DRMs and location. Figure 11. View largeDownload slide A νFν comparison of SEDs used for the integrated spectrum of GRB 080916C. The Band and Band+BB spectra are very similar due to the fact that the blackbody component contributes little to the spectrum. This indicates that the blackbody's contribution in Guiriec et al. (2015) may be due to neglecting the uncertainty of the DRMs and location. Table 4. GRB 080916C time-integrated balrog parameters (95 per cent HPD credible interval). Model  RA  Dec.  Ep  α  β  kT  Break scale  log  Z  Banda  132.8 : 136.5  −51.3 : −48.7  276.6 : 370.2  −0.91 : −0.82  −2.24 : −1.94  …  …  −2831.4  Band  120.9 : 144.8  −59.5 : −36.0  296.4 : 416.7  −0.95 : −0.86  −2.28 : −1.94  …  …  −1881.4  Band+BB  119.7 : 145.1  −60.2 : −36.0  238.6 : 422.9  −0.96 : −0.74  −2.26 : −1.92  1.0 : 325.0  …  −1881.1  SBPL  114.8 : 144.3  −62.2 : −37.0  136.2 : 194.3  −1.07 : −1.00  −2.07 : −1.82  …  1.78E-03 : 3.51E-01  −1878.8  Model  RA  Dec.  Ep  α  β  kT  Break scale  log  Z  Banda  132.8 : 136.5  −51.3 : −48.7  276.6 : 370.2  −0.91 : −0.82  −2.24 : −1.94  …  …  −2831.4  Band  120.9 : 144.8  −59.5 : −36.0  296.4 : 416.7  −0.95 : −0.86  −2.28 : −1.94  …  …  −1881.4  Band+BB  119.7 : 145.1  −60.2 : −36.0  238.6 : 422.9  −0.96 : −0.74  −2.26 : −1.92  1.0 : 325.0  …  −1881.1  SBPL  114.8 : 144.3  −62.2 : −37.0  136.2 : 194.3  −1.07 : −1.00  −2.07 : −1.82  …  1.78E-03 : 3.51E-01  −1878.8  Note. aFit including NaI 6 and 7. View Large We now fit the same data assuming no uncertainty in location, i.e. with a fixed DRM corresponding to the known location. Table 5 shows the recovered parameters corresponding to the fixed DRM fits. There are two important observations: the fitted parameters change when using a fixed DRM resulting in softer low-energy slopes and higher Ep and the Band+BB fit results in the well-known kT ∼ 40 keV blackbody with preferential marginal likelihood (Δlog Z ≃ 6) for the model over the Band function. This demonstrates that an overly precise location can induce apparent extra spectral components to compensate. Table 5. GRB 080916C time-integrated fixed DRM parameters (95 per cent HPD credible interval). Model  Ep  α  β  kT  Break scale  log  Z  Band  413.2 : 548.4  −1.02 : −0.95  −2.29 : −2.00  …  …  −237.4  Band+BB  955.9 : 1545.8  −1.28 : −1.16  −2.89: −2.06  41.2 : 48.64  …  −231.2  SBPL  205.1: 249.3  −1.09 : −1.11  −2.13 : −1.88  …  5.21E-03 : 3.98E-01  −234.4  Model  Ep  α  β  kT  Break scale  log  Z  Band  413.2 : 548.4  −1.02 : −0.95  −2.29 : −2.00  …  …  −237.4  Band+BB  955.9 : 1545.8  −1.28 : −1.16  −2.89: −2.06  41.2 : 48.64  …  −231.2  SBPL  205.1: 249.3  −1.09 : −1.11  −2.13 : −1.88  …  5.21E-03 : 3.98E-01  −234.4  View Large 5.2 Time-resolved analysis The high intensity of GRB 080916C caused Fermi to repoint during the burst (Abdo et al. 2009). The balrog assumed one spacecraft orientation when fitting spectra and locations, though we have designed it with the ability to fit multiple spectra simultaneously all sharing the same GRB location. We do not test this capability here and leave it for a future study. However, the 71 s duration used in the time-integrated analysis may introduce systematics as the spacecraft is repointing and a single spacecraft orientation is not valid. We can instead look at time-resolved spectra and use the spacecraft orientation corresponding to individual intervals to analyse the spectra. Ideally, all intervals should give similar locations. Using the same detector selection as was used in the time-integrated analysis, we choose the first eight intervals identified in the intriguing fine-time spectroscopy of Guiriec et al. (2015) and fit them with Band and CPL+BB SEDs. We additionally fit with the new GRB model identified in Guiriec et al. (2015) that consists of three spectral components, a CPL, a blackbody and a power law, where the CPL and PL have their spectral indices fixed to −0.7 and −1.5, respectively. The results from each time interval are detailed in Table 6. As with the time-integrated analysis, there was no significant Bayesian evidence for extra components in the spectra and a Band function fit the spectra sufficiently. The first interval provided an incorrect location as NaI 0 was not pointed close enough to the true location until after Fermi repointed. Therefore, we did not use NaI 0 in the first interval. Interestingly, while the marginal likelihood did not support the three-component model of Guiriec et al. (2015), the model resulted in inaccurate locations in several time intervals (see Fig. 12). Because of its fixed low-energy slope, the model's flexibility to fit the data is passed to the freedom in the location parameters, which adjusts to the poor fitting of the spectral model by giving an incorrect location. The marginal distributions of the Band fit and the CPL+BB+PL (Figs 13 and 14) demonstrate that while for the Band function properly constrained distributions are obtained, the more complex CPL+BB+PL model's parameters fill out their priors and hence lead to lower Bayesian evidence. One could argue that the incredibly small time intervals do not allow for more complex models to be statistically significant, however, the inability of the model to locate the GRB also makes it difficult to explain it as the true spectrum. However, recent work in Guiriec et al. (2016) shows that the model does help to explain the optical data of some GRBs, and hence, further investigation is required. This model is not able to explain the data of GRB 080916C. Figure 12. View largeDownload slide Locations produced from fitting the CPL+BB+PL function (left) and the Band model (right) to the interval T0+3.8–4.3 of GRB 080916C. The CPL+BB+PL model provides an inaccurate location. Figure 12. View largeDownload slide Locations produced from fitting the CPL+BB+PL function (left) and the Band model (right) to the interval T0+3.8–4.3 of GRB 080916C. The CPL+BB+PL model provides an inaccurate location. Figure 13. View largeDownload slide The marginal distribution for the Band (blue) and CPL+BB+PL (red) fit for the time interval T0+3.8–4.3 for GRB 080916C. K represents normalizations for the various spectral components. It is evident that for the CPL+BB+PL model, the CPL component is consistent with zero, allowing the BB to dominate the spectral peak. The differences in the location parameters are also clearly demonstrated. Figure 13. View largeDownload slide The marginal distribution for the Band (blue) and CPL+BB+PL (red) fit for the time interval T0+3.8–4.3 for GRB 080916C. K represents normalizations for the various spectral components. It is evident that for the CPL+BB+PL model, the CPL component is consistent with zero, allowing the BB to dominate the spectral peak. The differences in the location parameters are also clearly demonstrated. Figure 14. View largeDownload slide A νFν comparison between the Band and CPL+BB+PL for the interval T0+3.8–4.3. The width of the curves demonstrates the 95 per cent HPD marginal distributions including the error associated with location and, implicitly, the calibration. Figure 14. View largeDownload slide A νFν comparison between the Band and CPL+BB+PL for the interval T0+3.8–4.3. The width of the curves demonstrates the 95 per cent HPD marginal distributions including the error associated with location and, implicitly, the calibration. Table 6. GRB 080916C time-resolved parameters (95 per cent HPD credible interval). Time from T0  Model  Ep (keV)  α  β  kT (keV)  δ  RA  Dec.  log Z  −0.1 : 0.7  Band  154.6 : 618.9  −0.97 : −0.41  −4.86 : −1.89  …  …  82.5 : 157.9  −76.1 : −9.9  −346.0    CPL+BB  169.5 : 643.2  −0.97 : −0.43  …  1.00 : 314.79  …  89.1 : 158.6  −76.2 : −12.5  −346.4    CPL+BB+PL  22.2 : 97694.0  −0.7 (fix)  …  44.94 : 65.33  −1.5 (fix)  78.6 : 160.0  −75.1 : −7.1  −348.8  0.7 : 1.2  Band  169.5 : 425.5  −0.76 : −0.32  −5.00 : −2.18  …  …  81.4 : 154.9  −71.3 : −18.8  −326.2    CPL+BB  183.9 : 419.7  −0.76 : −0.33  …  1.00 : 217.38  …  81.6 : 155.6  −71.4 : −17.5  −326.5    CPL+BB+PL  16.4 : 82608.5  −0.7 (fix)  …  49.67 : 66.65  −1.5 (fix)  75.7 : 158.9  −71.3 : −6.2  −332.6  1.2 : 1.7  Band  625.8 : 1710.3  −0.97 : −0.71  −5.00 : −2.41  …  …  117.7 : 158.6  −64.7 : −14.3  −348.4    CPL+BB  169.9 : 2347.1  −0.99 : −0.40  …  6.77 : 797.42  …  110.0 : 158.2  −68.1 : −16.1  −347.6    CPL+BB+PL  17.4 : 98087.7  −0.7 (fix)  …  47.23 : 78.55  −1.5 (fix)  156.7 : 164.3  −5.7 : 9.2  −364.7  1.7 : 2.3  Band  93.2 : 255.7  −0.64 : −0.02  −4.20 : −1.77  …  …  92.1 : 155.9  −67.7 : −16.4  −399.3    CPL+BB  119.0 : 378.0  −0.71 : −0.19  …  1.01 : 569.82  …  96.1 : 154.7  −67.1 : −17.8  −400.1    CPL+BB+PL  15.0 : 88705.0  −0.7 (fix)  …  40.66 : 52.89  −1.5 (fix)  86.8 : 160.4  −66.6 : −0.8  −407.0  2.3 : 2.8  Band  270.2 : 634.1  −0.89 : −0.58  −5.00 : −2.45  …  …  100.2 : 154.5  −70.4 : −26.7  −336.0    CPL+BB  126.3 : 1513.9  −1.07 : −0.40  …  1.25 : 373.46  …  105.1 : 154.6  −69.5 : −26.5  −335.6    CPL+BB+PL  10.1 : 53433.7  −0.7 (fix)  …  47.46 : 60.70  −1.5 (fix)  74.7 : 153.8  −72.9 : −20.8  −346.8  2.8 : 3.4  Band  238.9 : 549.0  −0.86 : −0.54  −5.00 : −2.12  …  …  110.3 : 154.9  −66.2 : −21.2  −396.5    CPL+BB  267.9 : 566.2  −0.86 : −0.56  …  1.00 : 263.49  …  110.3 : 154.3  −65.1 : −20.4  −396.8    CPL+BB+PL  11.8 : 67813.5  −0.7 (fix)  …  52.77 : 68.14  −1.5 (fix)  91.7 : 158.5  −67.8 : −9.5  −401.8  3.4 : 3.8  Band  78.8 : 419.6  −0.87 : 0.02  −4.68 : −1.71  …  …  103.5 : 163.0  −66.2 : 2.2  −256.2    CPL+BB  86.2 : 767.5  −0.96 : −0.25  …  1.00 : 353.30  …  113.5 : 163.0  −63.0 : 0.4  −256.6    CPL+BB+PL  12.3 : 62089.7  −0.7 (fix)  …  36.52 : 55.27  −1.5 (fix)  151.8 : 163.0  −11.1 : 8.8  −262.5  3.8 : 4.3  Band  248.5 : 871.2  −0.94 : −0.57  −4.60 : −1.81  …  …  88.0 : 154.3  −68.6 : −20.0  −338.7    CPL+BB  227.9 : 2254.0  −1.02 : −0.57  …  1.22 : 808.97  …  93.1 : 153.9  −68.5 : −21.7  −339.1    CPL+BB+PL  10.5 : 57624.8  −0.7 (fix)  …  47.41 : 75.23  −1.5 (fix)  78.2 : 161.7  −68.0 : 3.6  −348.1  Time from T0  Model  Ep (keV)  α  β  kT (keV)  δ  RA  Dec.  log Z  −0.1 : 0.7  Band  154.6 : 618.9  −0.97 : −0.41  −4.86 : −1.89  …  …  82.5 : 157.9  −76.1 : −9.9  −346.0    CPL+BB  169.5 : 643.2  −0.97 : −0.43  …  1.00 : 314.79  …  89.1 : 158.6  −76.2 : −12.5  −346.4    CPL+BB+PL  22.2 : 97694.0  −0.7 (fix)  …  44.94 : 65.33  −1.5 (fix)  78.6 : 160.0  −75.1 : −7.1  −348.8  0.7 : 1.2  Band  169.5 : 425.5  −0.76 : −0.32  −5.00 : −2.18  …  …  81.4 : 154.9  −71.3 : −18.8  −326.2    CPL+BB  183.9 : 419.7  −0.76 : −0.33  …  1.00 : 217.38  …  81.6 : 155.6  −71.4 : −17.5  −326.5    CPL+BB+PL  16.4 : 82608.5  −0.7 (fix)  …  49.67 : 66.65  −1.5 (fix)  75.7 : 158.9  −71.3 : −6.2  −332.6  1.2 : 1.7  Band  625.8 : 1710.3  −0.97 : −0.71  −5.00 : −2.41  …  …  117.7 : 158.6  −64.7 : −14.3  −348.4    CPL+BB  169.9 : 2347.1  −0.99 : −0.40  …  6.77 : 797.42  …  110.0 : 158.2  −68.1 : −16.1  −347.6    CPL+BB+PL  17.4 : 98087.7  −0.7 (fix)  …  47.23 : 78.55  −1.5 (fix)  156.7 : 164.3  −5.7 : 9.2  −364.7  1.7 : 2.3  Band  93.2 : 255.7  −0.64 : −0.02  −4.20 : −1.77  …  …  92.1 : 155.9  −67.7 : −16.4  −399.3    CPL+BB  119.0 : 378.0  −0.71 : −0.19  …  1.01 : 569.82  …  96.1 : 154.7  −67.1 : −17.8  −400.1    CPL+BB+PL  15.0 : 88705.0  −0.7 (fix)  …  40.66 : 52.89  −1.5 (fix)  86.8 : 160.4  −66.6 : −0.8  −407.0  2.3 : 2.8  Band  270.2 : 634.1  −0.89 : −0.58  −5.00 : −2.45  …  …  100.2 : 154.5  −70.4 : −26.7  −336.0    CPL+BB  126.3 : 1513.9  −1.07 : −0.40  …  1.25 : 373.46  …  105.1 : 154.6  −69.5 : −26.5  −335.6    CPL+BB+PL  10.1 : 53433.7  −0.7 (fix)  …  47.46 : 60.70  −1.5 (fix)  74.7 : 153.8  −72.9 : −20.8  −346.8  2.8 : 3.4  Band  238.9 : 549.0  −0.86 : −0.54  −5.00 : −2.12  …  …  110.3 : 154.9  −66.2 : −21.2  −396.5    CPL+BB  267.9 : 566.2  −0.86 : −0.56  …  1.00 : 263.49  …  110.3 : 154.3  −65.1 : −20.4  −396.8    CPL+BB+PL  11.8 : 67813.5  −0.7 (fix)  …  52.77 : 68.14  −1.5 (fix)  91.7 : 158.5  −67.8 : −9.5  −401.8  3.4 : 3.8  Band  78.8 : 419.6  −0.87 : 0.02  −4.68 : −1.71  …  …  103.5 : 163.0  −66.2 : 2.2  −256.2    CPL+BB  86.2 : 767.5  −0.96 : −0.25  …  1.00 : 353.30  …  113.5 : 163.0  −63.0 : 0.4  −256.6    CPL+BB+PL  12.3 : 62089.7  −0.7 (fix)  …  36.52 : 55.27  −1.5 (fix)  151.8 : 163.0  −11.1 : 8.8  −262.5  3.8 : 4.3  Band  248.5 : 871.2  −0.94 : −0.57  −4.60 : −1.81  …  …  88.0 : 154.3  −68.6 : −20.0  −338.7    CPL+BB  227.9 : 2254.0  −1.02 : −0.57  …  1.22 : 808.97  …  93.1 : 153.9  −68.5 : −21.7  −339.1    CPL+BB+PL  10.5 : 57624.8  −0.7 (fix)  …  47.41 : 75.23  −1.5 (fix)  78.2 : 161.7  −68.0 : 3.6  −348.1  View Large This raises an interesting paradigm for model selection in GRBs, when the true location is known externally, models, such as the new physical models being fit to data (e.g. Burgess et al. 2014; Vurm & Beloborodov 2016), can be rejected if they do not provide proper location parameters. In future work we will address this with the synchrotron model of Burgess et al. (2014). 6 COMMENT ON THE CASE OF GW  150914 GBM The purported detection of electromagnetic radiation by GBM of Laser Interferometer Gravitational-Wave Observatory (LIGO) event GW  150914 (Abbott et al. 2016a,b; Connaughton et al. 2016) was in contrast with upper limits from the INTEGRAL-SPI/Anti-Coincidence Shield (ACS; Savchenko et al. 2016) and an independent analysis of the GBM data (Greiner et al. 2016). Our analysis in Greiner et al. (2016), which concluded that the data were consistent with purely background, differs from that described in Connaughton et al. (2016) in two ways: we used a different statistical model; and, possibly more importantly, we used only two detectors (NaI 5 and BGO 0) rather than all 14. We have found in our analysis herein that the use of all NaI detectors can result in systematically off locations owing to possible deficiencies in the GBM DRMs. Here, we use the balrog with two data sets: NaI 5 and BGO 0; and all 12 NaI and two BGO detectors. We use PL and CPL SED models to match Connaughton et al. (2016). Fig. 15 shows the results of the analysis when using all detectors with both SEDs. Both 1σ and 2σ regions of the GBM location intersect the LIGO arc. Disregarding the fact that using all 14 GBM detectors has the potential to introduce systematics, we can compare the fluences derived from these fits with the upper limits derived in Savchenko et al. (2016). Marginalizing over the location, Fig. 16 shows the marginal fluence distributions over the interval 70–2000 keV for both spectral functions. The PL function does not have a high-energy cut-off and therefore has less believable limits, although the 95 per cent HPD credible region does include the ACS upper limit. The CPL does have a high-energy cut-off and is clearly violating the ACS upper limits. The fact that systematics exist when using all detectors could propagate into the method of Blackburn et al. (2015) that relies on model rates in all GBM detectors built from the DOL templates discussed herein to estimate the significance of transient events in untriggered GBM data. If such systematics exist, and the event has a different spectrum than assumed from the spectral templates, then it is possible that the significance estimated for the GBM event is incorrect. Figure 15. View largeDownload slide The location of the event associated with GW  150914 using all GBM detectors and a PL photon (left) and CPL SED (right). The black points represent positions along the LIGO arc. Using all detectors results in location contours similar to Connaughton et al. (2016), yet it is possible that these are plagued with systematics. Figure 15. View largeDownload slide The location of the event associated with GW  150914 using all GBM detectors and a PL photon (left) and CPL SED (right). The black points represent positions along the LIGO arc. Using all detectors results in location contours similar to Connaughton et al. (2016), yet it is possible that these are plagued with systematics. Figure 16. View largeDownload slide The marginalized fluence distributions resulting from fits using all GBM detectors to the event associated with GW  150914. The red lines indicate the 95 per cent HPD credible interval for each distribution. The distributions are at odds with the ACS upper limits (Savchenko et al. 2016). Figure 16. View largeDownload slide The marginalized fluence distributions resulting from fits using all GBM detectors to the event associated with GW  150914. The red lines indicate the 95 per cent HPD credible interval for each distribution. The distributions are at odds with the ACS upper limits (Savchenko et al. 2016). Fig. 17 shows the results of fits obtained using data from the NaI 5 and BGO 0 detectors. The location is clearly unconstrained. Additionally, the spectrum is unconstrained. These results combined with the non-detection results of Greiner et al. (2016) point to a background (non-astrophysical) origin of the GBM counts data at the time of the GW  150914. In Greiner et al. (2016), we used two different methods to access the spectra with NaI 5 and BGO 0. We used a Poisson–Gaussian likelihood for fixed positions with maximum likelihood estimation and a Bayesian method where the background and source are fit simultaneously. Both methods found that the signal was consistent with background, the method presented here is similar to the maximum likelihood method of Greiner et al. (2016) except that the location parameters are free. Figure 17. View largeDownload slide The location of the event associated with GW  150914 using only NaI 5 and BGO 0 and a PL photon (left) and CPL SED (right). The black points represent positions along the LIGO arc. This subset of detectors has the best chance of representing the true location as they are those with the best viewing angle to the purported source (Connaughton et al. 2016). Figure 17. View largeDownload slide The location of the event associated with GW  150914 using only NaI 5 and BGO 0 and a PL photon (left) and CPL SED (right). The black points represent positions along the LIGO arc. This subset of detectors has the best chance of representing the true location as they are those with the best viewing angle to the purported source (Connaughton et al. 2016). In light of a new study by Veres et al. (2016), where a CPL was fit to the 14 detector GBM data to understand the plausible physical GRB scenarios that could explain the signal, we re-examine the spectral models of the GBM event. The authors use a method, which is not described, to incorporate multiple locations along the LIGO arc and the Monte Carlo maximum likelihood profile of the parameters. In that study, two sets of these Monte Carlos were created. Each set holds either the νFν peak constant as 1 MeV or the low-energy index at α = −0.42. We repeat this analysis to examine the marginal distributions under the balrog scheme. The locations of each fit are consistent with the LIGO arc at the 2σ level (see Fig. 18). We arrive at slightly different marginal distributions (Fig. 19) most likely due to the fact that we fully marginalize over all locations, though it is difficult to compare the methods without a detailed description from Veres et al. (2016). However, both fits violate the ACS upper limits shown in Fig. 20. Figure 18. View largeDownload slide The location of the event associated with GW  150914 using all GBM detectors repeating the analysis of Veres et al. (2016) where we hold Ep = 1 MeV fixed (left) and α = −0.42 fixed (right). The locations are marginally consistent with the LIGO arc at the 2σ level. Figure 18. View largeDownload slide The location of the event associated with GW  150914 using all GBM detectors repeating the analysis of Veres et al. (2016) where we hold Ep = 1 MeV fixed (left) and α = −0.42 fixed (right). The locations are marginally consistent with the LIGO arc at the 2σ level. Figure 19. View largeDownload slide The marginal distributions of the fixed low-energy index (purple) and fixed Ep (green) from the repeated analysis of Veres et al. (2016). The dotted lines represent the fixed values of each fit that show that the two different fits are inconsistent with each other. The fluence (erg cm−2) is integrated over 10–1000 keV. Figure 19. View largeDownload slide The marginal distributions of the fixed low-energy index (purple) and fixed Ep (green) from the repeated analysis of Veres et al. (2016). The dotted lines represent the fixed values of each fit that show that the two different fits are inconsistent with each other. The fluence (erg cm−2) is integrated over 10–1000 keV. Figure 20. View largeDownload slide The fluences from the fits of the CPL models following the analysis of Veres et al. (2016) compared with the ACS upper limit derived in Savchenko et al. (2016). Figure 20. View largeDownload slide The fluences from the fits of the CPL models following the analysis of Veres et al. (2016) compared with the ACS upper limit derived in Savchenko et al. (2016). 7 DISCUSSION Localizing GRBs is the most important function of GBM in the era of multimessenger astronomy. With the possibility of joint LIGO–GBM observations in the near future, GBM localizations have the potential to drastically shrink the broad localization arcs provided by LIGO leading to higher probability that imaging instruments can zero in on the location of a LIGO trigger. With this goal in mind, we have proposed a modern method to localize GBM GRBs that simultaneously fits the spectrum and location of GRBs via Bayesian nested sampling. The technique eliminates the systematics inherent in the standard GBM localization technique reported in Connaughton et al. (2015), though some issues still remain that may be the result of poor spacecraft modelling or improper spectral model choice. Proper spacecraft modelling can be addressed with better simulations but model choice is a problem that plagues GRB spectroscopy for even the very wise cannot see all ends. In attempting to find the true location of GRBs with systematic offsets via balrog, we found that special care must be taken with both model and detector selection. The following methodology was developed. A subset of detectors with 5σ significance above background (Li & Ma 1983) and less than 60° angular separation from the GBM flight software position are initially used and a location is determined using the CPL model. Once a candidate location has been determined, it is checked that the detectors selected view the final location. If not, the selection is modified and repeated. This process is repeated while also iterating model choice until the highest marginal likelihood is found that we then take to correspond to the true location.Proper location depends heavily on detector selection is indicative of deficiencies in the detector response. This could be explained if: (i) the >60° off-axis response of the individual detectors has some calibration uncertainties; (ii) the Bayesian fitting scheme does not perfectly cover the noise-only detectors; (iii) the simulations include deficiencies because they are based on imperfect mass models of the detectors and the spacecraft and the corresponding spacecraft-related scattering. The effect of these deficiencies has been noted already at the pre-launch source survey, when the full Fermi spacecraft with all GBM detectors mounted was irradiated with a radioactive source. The simulations were never able to reproduce the measurements, in particular for detectors that are at large angles relative to the source beam. Since these mismatches did not show up for the calibration measurements of individual detectors (Bissaldi et al. 2009), we presently suspect the above effect (iii) to be dominant. A deeper investigation into these deficiencies may alleviate the need for careful detector selection that currently makes using the balrog difficult in real-time follow-up. Note that this study only analyses a small fraction of the GRBs listed in Connaughton et al. (2015). However, these are the tail cases (those corresponding to the best and worst localizations made by the DOL) and balrog performs well in both. It is conceivable that for the median of the distribution, balrog will perform equally well. We advocate for a deeper study to access this, but feel that the introduction of the concept herein provides the required impetus to modernize GRB location algorithms. It is also worth restating that this method is a superset of the original DOL method. balrog is computationally expensive and fast location estimates will require its implementation on high-performance computing clusters. A standard location on TRIGDAT data takes ∼1–5 min on a 48-core node without hyperthreading. The main computational bottleneck is the DRM generation. In the future, it may be possible to pre-compute a finer base DRM grid, but the multiplicity of spacecraft orientations that require recomputing of the atmospheric scattering will make such a task difficult. A secondary finding of this study is the effect of including the location errors in spectra. The spectral shape of GRB 080916C for both time-integrated and time-resolved analysis favour a simple Band function rather than multiple components. The new model proposed in Guiriec et al. (2015) is not favoured and results in incorrect locations due to its fixed low-energy slope. To reject a spectral model in this way, it is required that a GRB is localized externally by an instrument such as Swift. In the future, optical/radio counterparts will be much easier to identify given the ∼few deg2 error regions as compared to the ∼100 deg2 error regions (where afterglow identification was possible already in about 50 per cent of the cases; Singer et al. 2013, 2015; Lipunov et al. 2016). While these findings focus on a single GRB, they provide impetus to study these effects on a broader sample that will be the subject of a future work. We reiterate that this opens a new path for GRB model selection. While the Band function's flexibility is likely sufficient to allow for proper location, the push for using physical models in GRB spectroscopy that have less flexibility will allow for quick rejection of these models if they result in inaccurate locations. The systematics identified when using all NaI detectors support the claims of Greiner et al. (2016) regarding the purported signal associated with GW  150914 (Connaughton et al. 2016). While including all NaI detectors results in a location that touches the LIGO arc, it also includes most of the Earth. When using a detector selection that corresponds to the position of the LIGO arc, the location contours encompass the entire sky. This, in combination with the lack of signal found in Greiner et al. (2016), places serious constraints on the astrophysical origin of the event. Furthermore, the balrog can be used for future LIGO events as it provides accurate locations for GBM transients that can be employed in optical follow-up strategies. With these findings, we advocate for adopting the balrog for the location of GBM GRBs as it provides realistic and accurate locations suitable for follow-up strategies. Additionally, the shifts in spectral properties found in GRB 080916C warrant an adoption of this scheme for spectral analysis of GBM GRBs but could also be used for similar instruments. We will investigate the GBM GRB catalogues in a forthcoming publication. Such a scheme could be implemented in modern analysis tools such as the Multi-Mission Maximum Likelihood framework (3ML; Vianello et al. 2015) that would allow the balrog to be used to fit data simultaneously with other instruments while still incorporating the uncertainty in location. ACKNOWLEDGEMENTS This work made use of several open source python packages and we are thankful to the authors of matplotlib (Hunter 2007), astropy (Astropy Collaboration 2013), seaborn (Waskom et al. 2016) and pymultinest (Buchner et al. 2014). The authors are thankful to the authors of the original GBM DRM generation code that helped in designing balrog and to Gudlaugur Jóhannesson for insight into making the code more efficient. Additional statistical discussion with Brandon Anderson and Giacomo Vianello helped steer this work in the right direction. We thank Andreas von Kienlin (MPE) and Marc Kippen (Los Alamos) for discussion on the GBM response calibration and simulation. Footnotes 1 GBM and BATSE DRMs differ from classical instrument DRMs because the effective area and response are combined into one matrix. 2 Our balrog tool is completely unrelated, except in name, to the BaLROG project to examine the dynamical influence of bars in galaxies described by Seidel et al. (2015). 3 The actual expressions given in the pgstat section of Arnaud et al. (2015) are −2 times the logarithm of the likelihoods given here. There are also some terms with no S dependence that have been omitted here and the ambiguous ± sign in equation (10) has been made definite. REFERENCES Abbott B. P. et al.  , 2016a, Phys. Rev. Lett. , 116, 061102 CrossRef Search ADS   Abbott B. P. et al.  , 2016b, ApJ , 826, L13 CrossRef Search ADS   Abdo A. A. et al.  , 2009, Science , 323, 1688 CrossRef Search ADS PubMed  Arnaud K., Dorman B., Gordon C., 2015, Xspec User's Guide, v12.9.0, available at: https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html Astropy Collaboration, 2013, A&A , 558, A33 CrossRef Search ADS   Band D., Matteson J., Ford L., 1993, ApJ , 413, 281 CrossRef Search ADS   Bissaldi E. et al.  , 2009, Exp. Astron. , 24, 47 CrossRef Search ADS   Blackburn L., Briggs M. S., Camp J., Christensen N., Connaughton V., Jenke P., Remillard R. A., Veitch J., 2015, ApJS , 217, 8 CrossRef Search ADS   Briggs M. S. et al.  , 1996, ApJ , 459, 40 CrossRef Search ADS   Buchner J. et al.  , 2014, A&A , 564, A125 CrossRef Search ADS   Burgess J. M., 2014, MNRAS , 445, 2589 CrossRef Search ADS   Burgess J. M. et al.  , 2014, ApJ , 784, 17 CrossRef Search ADS   Connaughton V. et al.  , 2015, ApJS , 216, 32 CrossRef Search ADS   Connaughton V. et al.  , 2016, ApJ , 826, L6 CrossRef Search ADS   Feroz F., Hobson M. P., Bridges M., 2009, MNRAS , 398, 1601 CrossRef Search ADS   Gelman A., Carlin J. B., Stern H. S., Dunson D. B., Vehtari A., Rubin D. B., 2014, Bayesian Data Analysis, 3rd edn. CRC Press, Boca Raton, FL Goldstein A. et al.  , 2012, ApJS , 199, 19 CrossRef Search ADS   Greiner J., Burgess J. M., Savchenko V., Yu H. F., 2016, ApJ , 827, L38 CrossRef Search ADS   Guiriec S. et al.  , 2015, ApJ , 807, 148 CrossRef Search ADS   Guiriec S. et al.  , 2016, ApJ , 831, L8 CrossRef Search ADS   Hunter J. D., 2007, Comput. Sci. Eng. , 9, 90 CrossRef Search ADS   Kippen R. M. et al.  , 2007, in Ritz S., Michelson P., Meegan C. A., eds, AIP Conf. Proc. Vol. 921, The First Glast Symposium. Am. Inst. Phys., New York, p. 590 Li T. P., Ma Y. Q., 1983, ApJ , 272, 317 CrossRef Search ADS   Lipunov V. M. et al.  , 2016, MNRAS , 455, 712 CrossRef Search ADS   Loredo T. J., Lamb D. Q., 1992, in Paciesas W. S., Fishman G. J., eds, AIP Conf. Proc. Vol. 265, Gamma-Ray Bursts. Am. Inst. Phys., New York, p. 414 Meegan C. et al.  , 2009, ApJ , 702, 791 CrossRef Search ADS   Paciesas W. S., Pendelton G. N., Wilson R. B., Fishman G. J., Meegan C. A., 1988, BAAS, 726 Savchenko V. et al.  , 2016, ApJ , 820, L36 CrossRef Search ADS   Seidel M. K., Falcón-Barroso J., Martínez-Valpuesta I., Díaz-García S., Laurikainen E., Salo H., Knapen J. H., 2015, MNRAS , 451, 936 CrossRef Search ADS   Singer L. P. et al.  , 2013, ApJ , 776, L34 CrossRef Search ADS   Singer L. P. et al.  , 2015, ApJ , 806, 52 CrossRef Search ADS   Skilling J., 2004, in Fischer R., Preuss R., Toussaint U. V., eds, AIP Conf. Proc. Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering. Am. Inst. Phys., New York, p. 395 van Dyk D. A., Connors A., Kashyap V. L., 2001, ApJ , 548, 224 CrossRef Search ADS   Veres P., Preece R. D., Goldstein A., Mészáros P., Burns E., Connaughton V., 2016, ApJ , 827, L34 CrossRef Search ADS   Vianello G. et al.  , 2015, preprint (arXiv:1507.08343) Vurm I., Beloborodov A. M., 2016, ApJ , 831, 175 CrossRef Search ADS   Waskom M. et al.  , 2016, seaborn: v0.7.0 (January 2016) , available at: http://dx.doi.org/10.5281/zenodo.45133 APPENDIX A: LOCATION PLOTS Figure A1. View largeDownload slide The plots of the remaining nine GRBs well located by the DOL and balrog. Figure A1. View largeDownload slide The plots of the remaining nine GRBs well located by the DOL and balrog. Figure A2. View largeDownload slide The plots of the remaining nine GRBs with large DOL systematic error. Figure A2. View largeDownload slide The plots of the remaining nine GRBs with large DOL systematic error. © 2017 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial