The EGRET Likelihood Analysis Program, LIKE

J.R. Mattox and D.J. Macomb

10-September-1993


TABLE OF CONTENTS

  1. Introduction
  2. The Likelihood of EGRET Data
  3. The Gamma-ray Model
  4. Maximum Likelihood Parameter Estimation
  5. Likelihood Ratio Test
    1. Test of EGRET Source Existence
    2. The Uncertainty of the Estimate
      1. The Uncertainty of the Counts Estimate
      2. Determination of a Counts Upper Limit
      3. The Uncertainty of the Point Source Location
  6. Using LIKE
    1. Adjusting the Model Parameters
    2. Evaluating the Model
    3. Likelihood Testing
    4. LIKE Output
  7. References
  8. Appendix A. A Sample LIKE Session
  9. Appendix B. Commands
  10. Appendix C. The CTL file

1  Introduction

The EGRET point source analysis system (LIKE) is used to examine the EGRET data for point gamma-ray sources. Specifically, it is designed to find an excess above the predicted diffuse gamma-ray flux which is consistent with the instrument point spread function (PSF). Two approaches have been developed for point source analysis: cross-correlation (Mayer-Hasselwander and Simpson 1988), and maximum likelihood. The cross-correlation method relies upon correlating the instrument PSF with the measured data and comparing this with the correlation between the PSF and the expected diffuse flux. LIKE uses the maximum likelihood approach (Cash, 1979) and has advantages over cross-correlation, especially in parameter estimation. This document is intended to describe the basics of the maximum likelihood approach as well as serve as a users guide to the functions of the LIKE program.

Experience with flight data indicates that the method of maximum likelihood provides reliable gamma-ray results. See Mattox et al. (1994) for a full description for the use of likelihood for EGRET data analysis. Monte Carlo characterization of the method is in progress (Willis et al., 1993). The method allows simultaneous estimate of background model parameters (a multiplier for a diffuse galactic emission model and the intensity of isotropic emission) and the intensity of emission from a hypothetical point source distributed as the EGRET PSF.

The maximum likelihood ratio test provides a means to evaluate the strength of the evidence for the existence of a point source of gamma rays. By Wilks' theorem, -2 times the log of the maximum likelihood ratio (the likelihood of the data without a point source, or the null hypothesis, divided by the likelihood of the data with a point source which is the alternative hypothesis) is distributed as c2 with one degree of freedom (source strength at a predetermined location). If the source is not detected, the dependence of the likelihood ratio on the intensity of the hypothetical point source may be used to obtain an upper limit. The dependence of the likelihood ratio on the assumed position of the point source allows for estimation of the source position and the uncertainty of this position. In addition, estimates of the intensity of emission for sub-intervals of the EGRET energy range are used to obtain EGRET spectra.

LIKE is a Fortran program which accomplishes all the above tasks. LIKE uses the counts and exposure maps produced by the EGRET programs MAPGEN and INTMAP respectively. In addition, it uses a map of expected diffuse galactic emission (e.g., the output of the EGRET team programs GALDIF and ISODIF). All input maps are expected to be in standard FITS format. Reading and writing of FITS files uses FITSIO routines. Section 6 of this document describes some of the basic capabilities of LIKE. Before this, however, sections 2-5 will describe the theoretical basis for the LIKE algorithms as well as some technical aspects of LIKE calculations.

2  The Likelihood of EGRET Data

An appropriate starting point is to describe the likelihood which is the joint probability of the observed EGRET data, for a specific model which will be denoted as q. The Poisson probability for each pixel in the binned gamma-ray counts map is
pij = qijnije-qij
nij!
where i and j index the pixels of the two-dimensional map, nij is the number of counts in the pixel, and qij is the number of counts predicted by the model. The likelihood is then the product of the pixel probabilities within a user-defined radius, RANAL - for Radius of ANALysis, of the hypothetical point source
Lq =
Õ
ij 
pij
The radius is typically chosen to be large enough to cover almost all the EGRET PSF so that all source counts are included as well as providing sufficient data to estimate the background model parameters, yet small enough to avoid increasing systematic errors involved in the estimated gamma-ray background model. The log of the likelihood is used in hypothesis testing and is usually more convenient calculationally. The log of the likelihood is then
lnLq = log(Lq) =
å
ij 
nijlog(qij)-
å
ij 
qij - å
log(nij!)
Since the last term is model independent, it is not relevant for the likelihood ratio and will not be considered further. The log of the likelihood thus reduces to
lnLq =
å
ij 
nijlog(qij)-
å
ij 
qij         (1)
which serves as the basis for LIKE hypothesis testing and parameter estimation based upon the measured counts and the gamma-ray model.

3  The Gamma-ray Model

The basic gamma-ray model, q, is constructed from three quantities; the counts from a possible point source, the diffuse galactic flux, and a constant diffuse extragalactic flux. The LIKE program will read a diffuse emission prediction map in FITS format (e.g. GALDIF and ISODIF output). The standard EGRET team prediction maps are based upon recent observations of galactic CO and HI and assumed cosmic ray interactions. It is worth noting that LIKE will also work without a diffuse prediction map, but the results become less and less reliable as the candidate source point approaches the galactic plane.

To convert the prediction map to counts, it is multiplied by the exposure for each pixel, Eij, to produce a predicted full resolution diffuse counts map G·ij. The LIKE program convolves this map with the appropriately calculated EGRET PSF to get the estimated diffuse galactic counts per bin

Gij =

å
kl 
G·kl PSF(Dijkl)


å
kl 
PSF(Dijkl)
where Dijkl is the distance between pixels indexed by ij and kl. Note that here, and elsewhere, we are ignoring the energy dependence of the PSF which has been integrated out in calculating the PSF for the measured energy range of the map in use. The summations are done for all kl pixels within the map which are less than RANAL from pixel ij and which have exposure.

With the addition of the term involving the estimated extragalactic diffuse flux multiplied by the exposure, the general pixel counts model becomes

qij = gm Gij+gb Eij +
å
k 
ck PSF(lk,bk,i,j)    (2)
where gm is Gmult, the galactic multiplier, and gb is Gbias, the isotropic diffuse flux (with units cm-2s-1 sr-1). The number of counts detected from the candidate source at lk,bk is ck (LIKE works equally well in celestial or galactic coordinates), and PSF(lk,bk,i,j) is the fraction of the PSF centered on lk,bk contained in pixel ij. Values of gm, gb, ck, lk,bk may be specified manually, or estimated by LIKE. This is the model used in the likelihood ratio test by testing the null hypothesis, ck = 0 (no source present), against the hypothesis that ck has the estimated value where gm and gb have their optimal values for both hypotheses. The estimation of the model parameters is described in the next section.

4  Maximum Likelihood Parameter Estimation

The estimate of a parameter is the value which maximizes the likelihood, or equivalently the log of the likelihood (eq.1). In the limit of a large number of counts, the estimated value of the parameter [^(r)] is normally distributed (Kendall & Stuart 1961, 18.14) with mean equal to the true value and standard deviation given by
s = é
ë
- 2lnL
r2
ù
û
-1/2
 
      (3)
The maximum of eq. 1 is obtained by LIKE for the optimization of ck with the DBRENT routine of Numerical Recipes (Press et al. 1986) which uses parabolic interpolation with information from the first derivative. For optimizing other parameters, or for simultaneous optimization of multiple parameters, Newton-Raphson iteration is used. This method is particularly efficient because lnL is a nearly parabolic function of ck, gm, and gb (which is the form assumed in the Newton-Raphson method). The use of the derivatives in analytical form drastically increases the efficiency of the calculations.

lnL
gm
=
å
ij 
Gij ( nij
qij
- 1)

lnL
gb
=
å
ij 
Eij ( nij
qij
- 1)

lnL
ck
=
å
ij 
PSF(lk,bk,i,j) ( nij
qij
- 1)

2 lnL
gm2
= -
å
ij 
nij( Gij
qij
)2

2 lnL
gb2
= -
å
ij 
nij( Eij
qij
)2

2 lnL
ck2
= -
å
ij 
nij( PSF(lk,bk,i,j)
qij
)2

2 lnL
gm gb
= 2 lnL
gb gm
= -
å
ij 
nij Gij Eij
qij2

2 lnL
ck gb
= 2 lnL
gb ck
= -
å
ij 
nij PSF(lk,bk,i,j) Eij
qij2

2 lnL
ck gm
= 2 lnL
gm ck
= -
å
ij 
nij PSF(lk,bk,i,j) Gij
qij2

When estimating gm, gb, or ck, LIKE reports a 1 s uncertainty using eq. 3. Except in the case of a counts upper limit calculation, LIKE constrains all parameters to be non-negative. LIKE may also be used to estimate a source position via the downhill simplex method (again Numerical Recipes, AMOEBA code). However, the derivative is not available as an analytical expression for determining source position errors. LIKE provides the option of making a fine map of the likelihood ratio test statistic (command MF or LE, described in section 6) to estimate the source position and uncertainty.

5  Likelihood Ratio Test

Once the model parameters have been evaluated, hypothesis testing can begin. The likelihood ratio test has proven to be a workable test with good properties for diverse purposes (Eadie et al. 1971; Kendall & Stuart). It was first applied to high-energy gamma-ray astronomy by Pollock et al. (1981 and 1985). The likelihood test statistic is
TS = 2( max  lnLq(q Î Q) - max  lnLq(q Î n))
where max  lnLq(q Î Q) is the maximum value of the log of the likelihood for the alternative hypothesis (a model within the parameter space Q). Conversely,

max  lnLq(q Î n) is the maximum value of the log of the likelihood for the null hypothesis (a model within the subspace n of Q; i.e. some parameters are restricted). The Likelihood Ratio Test statistic is especially useful because in the limit of a large number of counts (see Cash 1979 for an approximation of the distribution for small number events), the distribution of TS is c2l under the null hypothesis, where l is the difference in the number of free parameters between the alternative hypothesis, and the null hypothesis.

5.1  Test of EGRET Source Existence

A crucial application of the likelihood ratio test  is the examination of EGRET data to determine if a point source exists at a specific point. As stated above, the extra free parameter in the alternative hypothesis is the number of counts detected from the putative source, candidate source counts being fixed at zero for the null hypotheses. If the source counts are the only parameter varying, so gb and gm are optimized for both models and the source location is assumed to be known, the TS is distributed as c21 in the null hypothesis. Thus, the significance (pre-trials) of a point source is
S = ó
õ
¥

TS 
c21(x) dx
which gives the probability of the test statistic being at least as large as that calculated and may be found readily in graphical or tabular form (e.g. Particle Properties Data Booklet) or calculated with standard software (e.g. IDL, Numerical Recipes or Mathematica). A limited tabulation follows with the equivalent gaussian distribution ßigma". Figure 5.1 also plots the log of the significance for a limited range of the TS. The c2 statistic for one d.o.f. is nearly exponential as shown by the relatively straight line behavior of the log. It is worth noting that the double sided Gaussian, s, is equal to [ÖTS].

TSSignificancessingle sidedsdouble sided
10.320.481
2.0000.15731.0051.41
3.845.0×10-21.651.96
5.22.3×10-22.02.3
6.61.0×10-22.32.6
10.31.3×10-33.03.2
125×10-43.33.5
17.33.2×10-54.04.2
208×10-64.44.5
304×10-85.45.4
403×10-106.26.3
502×10-1277
1002×10-231010
50010-1102323
100010-2193232
300010-6545555

Table 5.1 Tabulation of Significance levels for some representative values of the TS.

The pre-trials log of the significance for a
restricted range of the TS.

Figure 5.1 The pre-trials log of the significance for a restricted range of the TS. larger
Theoretically, the two position coordinates for an unidentified source should also be considered to be degrees of freedom; necessitating a c2 statistic with 3 d.o.f. for testing source existence. In practice, however, it is found that using 1 d.o.f. and accounting for source location with a multiplicative factor (the number of trials) approximately equal to the EGRET beam width divided into the area of sky searched provides better results.

5.2  The Uncertainty of the Estimate

The likelihood ratio test  may be used to evaluate the uncertainty of the estimate of a model parameter. In this case, the equivalent null hypothesis has the parameter in question fixed at its true (but unknown) value. Since the distribution of the test statistic is c2l under the null hypothesis, the uncertainty of the parameter estimate [^(r)] which maximizes the likelihood may be readily obtained. With confidence (C = 1- S), S being the probability that the true value lies outside the specified range
C = 1- ó
õ
¥

h 
c2l(x)dx
the true parameter(s) reside in the region where \mid TS([^(r)])-TS(r)\mid £ h.

This may be applied to three important questions with LIKE: counts uncertainty, counts upper limit when the null hypothesis cannot be rejected, and position uncertainty determination.

5.2.1  The Uncertainty of the Counts Estimate

The confidence region is that for which TS ³ TS([^(r)])-h. A value of h = 1 corresponds to C = 68% (one sigma).

5.2.2  Determination of a Counts Upper Limit

If the candidate source is not detected, an upper limit can prove useful. This estimate is given by the number of counts ck for which TS(ck) = TS(ck = 0)-h is the confidence level C upper limit. If the unconstrained counts estimate is negative (which is non-physical), the upper limit is assumed to be the difference between the ``statistical'' upper limit and the counts estimate. (The ``statistical'' upper limit is the value of counts where lnL drops by 1.92 (for 95% confidence) from it's peak at the estimate.) This procedure is ad hoc and, while preliminary work shows good agreement, needs to be justified with further Monte Carlo work.

The confidence level versus eta for 
the chi-squared(2) distribution

Figure 5.2 The confidence level versus h for the c22 distribution. larger

5.2.3  The Uncertainty of the Point Source Location

In this case, the alternative hypothesis has two extra free parameters (i.e. lk,bk) while the source counts are fixed. Therefore, the TS under the null hypothesis is distributed as c22 in the limit of a large number of gamma-ray events (see Cash 1979 for an approximation of the distribution for small number events). The fine TS maps produced with LIKE MF command (see section 6.4) may be plotted with other software to show contour lines which correspond to 50,68,95, and 99% confidence regions (corresponding h, 1.4,2.3,6.0,9.1). Also, the LE command will make a plot of these contour lines and show the centroid and maximum TS positions. Figure 5.2 shows C vs h for l = 2(the c2 distribution for 2 d.o.f. is an exponential).

6  Using LIKE

The previous sections provided a summary of the concepts behind the maximum likelihood method on which the LIKE program is based. It is now time to concentrate on the specifics of running the program itself. LIKE runs as a stand alone program and is linked with FITSIO for much of its file i/o and with PGPLOT for graphics output. The user interface is command line driven and the output graphics can accommodate Xwindows, Postscript and other devices. The first step in running like is to have the EGRET data analysis environment variables setup correctly. If you have been using other EGRET software, most of this is probably done. If not, ask the system manager at the site running LIKE about what needs to be done to define your environment. The following list defines the environment variables which should be properly set for running LIKE.

This is a minimal set and some options may depend on other variables, but if LIKE does not find a file where it thinks it should be, the defaults can often be overridden by command line input. Properly setting these environment variables, however, will ensure that LIKE finds the proper counts, exposure, and calibration files.

The basic input to LIKE consists of standard EGRET counts and exposure maps and the EGRET calibration files. Maps up to 720 longitude and 360 latitude pixels can be used (enough for the entire sky at 0.5 degrees resolution). Note that LIKE does not do any summing of input mapfiles; co-adding maps should be done before running LIKE (e.g., using the EGRET team ADDMAP program). The standard Energy > 100 Mev maps can be retrieved by simply giving the four number viewing period (for example, period 5 Þ 0050). Similar shortcuts for other available maps (i.e. phase 1 total, phase 2 total, etc.) will exist. LIKE has no facility for perusing or adjusting the input maps so if one wants to print a FITS header or display a portion of the input map, use ancillary software such as the EGRET SHOW or SKYMAP programs.

The main calibration file required by like is the PSD file (Point Spread Function), although the SAR (Sensitive Area) and EDF (Energy Dispersion) files are used as well. Each calibration file contains information derived from pre-launch calibration at the Stanford Linear Accelerator Center (SLAC). For more information on EGRET calibration, refer to Thompson et. al., 1993. The PSD file contains the EGRET point spread function at 20 discrete energies, 3 azimuthal angles, and 9 inclinations. All calibration files can usually be found in the directory pointed to by the CALIB_DIR environment variable. There are several calibration files available which represent different operating conditions for the EGRET instrument, but the difference in the resulting PSF is not believed to be significant.

The actual EGRET PSF used by LIKE is constructed based upon the measured energy range of the map and will be a combination of the PSF's measured at discrete energies combined assuming a candidate source has some given power-law index, GAMMA. This also requires using the SAR and EDF calibration files since the sensitive area as a function of energy and the dispersion of measured gamma-ray energies are utilized. Averaging over inclination and azimuth are also performed. By default, the CTL file is searched for in the working directory with filename='CTL' and, if not found, a full pathname should be given. The CTL file itself sets up many LIKE parameters. The values of entries in the CTL file don't normally need to be adjusted, and if they are, should be adjusted with caution. An annotated CTL file is explained in Appendix 3.

Additional input to LIKE is the diffuse model map, necessary for the estimate of gmult, and the ancillary PSF file which is a list of EGRET sources which can be included in the analysis. LIKE will run without a diffuse model map. However, especially near the galactic plane, results will not be optimal. The input diffuse maps will note whether or not the map has already been convolved with the EGRET PSF for the appropriate energy range. The diffuse model map must match the selected energy range of the input counts map if it has been convolved with the PSF.

Once the setup stage is complete, LIKE is ready to run. The initial state of LIKE is that the Region of Interest (ROI), which is the active map area for several functions, is set to the entire map; RANAL, a critical parameter, which is the radius of analysis for determining model parameters has been input (the default is energy dependent), and g_b and g_m are set to their energy dependent initial values. The next step is to perform the analysis. Everyone using LIKE will have their own preferences and approaches to this process. There is no guide like experience! Even so, the following subsections describe a few of the major options of LIKE in a somewhat logical progression from setting up the parameters to outputting results (also see Appendix 1 for a sample LIKE session or Appendix 2 for a complete list of LIKE commands).

6.1  Adjusting the Model Parameters

A reasonable starting place is with the commands to alter the default parameters used in the analysis. These commands are likely to be used at the beginning as well as throughout your LIKE session. They fall into two main categories; changing the relevant map and changing some parameter of the current analysis.

Two of the most fundamental adjustments are to change the ROI or the energy range. These are accomplished with the S and E commands. Changing the energy (E)will allow you to select from a different range of a multi-range FITS file or to choose a completely new map; of course, a new PSF will also be calculated at the same time. This command essentially starts the LIKE session over again as all previously performed calculations are no longer valid for the new input map. The command to define a new ROI (S) (input the defining coordinates as lower LONG, upper LONG, lower LAT, upper LAT) is less drastic. The ROI is the portion of the input map which is the area of interest for many commands which search over large areas of the map (i.e., searching for a point source with the MS command described below) and also defines the mapping region for many output functions. Note that RANAL can extend out beyond the defined ROI with those points contributing to the analysis. After adjusting the ROI, many aspects of the current analysis may still be valid, but some calculated values may need to be changed. Alternate methods for setting the ROI are that a rectangle may be selected by specifying the rectangle center and width in either direction or that the coordinates of a point and the half-width of a square about the point may be input.

Adjustments to other parameters for the current analysis are usually accomplished through the Pnn command set. These commands are quite valuable and may effect the analysis in a major way. Perhaps the two most important are PA, which adjusts the ``other PSF map'', and PR which adjusts the range for the likelihood analysis (RANAL). The radius for which data is included in the likelihood test may be as large as 180 degrees. The PSF matrix only goes to 20 degrees - data beyond 20 degrees will pertain to the determination of Gmult and Gbias only for unmodeled sources. RANAL should be large enough to include most of the EGRET PSF (which is an energy dependent value; see Thompson et. al., 1993), yet small enough to keep errors in the background model from adversely effecting the analysis. As the size of RANAL is increased, however, more data is used for estimating the model parameters decreasing statistical error. This is an example of the classic tradeoff between systematic and statistical errors in data analysis. The correct value of RANAL has not been definitively prescribed.

In some cases, the field of view may contain more than one source. The analysis may be affected by wings of the PSF from any individual source affecting the likelihood estimates for a separate candidate source. This problem can be obviated by adjusting the model to include PSF's from any nearby, known sources with the PA command. Up to 500 PSF's from known sources can be accommodated. An input catalog of source positions and counts for sources may be used. The PA command invokes the PSF_ADJ sub-menu. The sub-command I will read such a disk file, and build the OTHER PSF component of the model from this list. At this point, each new PSF is active, or included in the model. The sub-command X will delete an individual PSF. The sub-command A will add the current active psf, while the sub-command D will add a modification of the active psf. The sub-command M will toggle the psf active Mark (for LM analysis). These commands will all prove valuable in adjusting the model for multiple sources in the field-of-view before likelihood testing or parameter estimation for an unidentified source.

Also notable is the PMS command, which toggles the SPECTRAL flag on and off, also changes the value of delta_lnL used in an upper limit evaluation to 1.0 (the appropriate value for 1 sigma error estimates - corresponding to a confidence of 84% - the integral of the gaussian distribution from -¥ to 1 s). Delta_lnL may be changed otherwise with the PML command. If the SPECTRAL flag is true, the counts obtained in the L or LM commands are written in a format which can be directly pasted into a SPECTRAL input file. See the command summary in Appendix 2 for information on other parameters which can be adjusted.

6.2  Evaluating the Model

Once LIKE is setup to the users satisfaction, the evaluation of model parameters for any position in the map can begin. This section briefly describes those LIKE commands which estimate parameters in the gamma-ray model (i.e., the source counts - Cxx commands, gas multiplier - Gxx commands, or gas bias -Bx commands). These commands are the most basic commands for estimating parameters.

Each parameter is estimated according to the methods described in Section 4. The user has the option of estimating each parameter individually, in pairs, or all 3 model parameters simultaneously. If a multi-parameter estimation fails to converge, the number of parameters being optimized is reduced. A warning is printed if only the reduced optimization is achieved. Turn on the verbose flag (PMV) to see exactly what transpires in such a case. Negative values of Counts, Gmult, & Gbias are not permitted. Two critical parameters are Gmult_nom, and Gbias_nom - the approximate values expected. These are used to start the Newton-Raphson method as close as possible to the optimal values in order to achieve reliable convergence. Both are set automatically upon loading a counts map, but may be changed with PMN. One can estimate any parameter with extant values of the other model parameters, or adjust the others manually. Counts units are actually an estimate of the number of photons received whereas gmult is a multiplier for the diffuse intensity model (nominally 1.0 for E > 100 MeV) and gbias is reported (or specified) in units of 10-5cm-2 s-1 sr-1 (Also nominally 1.0 for E > 100 MeV).

For example, estimation of source counts for a specific location with the extand values of gmult and gbias can be accomplished the CC command while the `C command alone allows the user to adjust gmult and gbias if desired. The CB command will evaluate counts and gbias simultaneously with a user set value of gmult while CBB will evaluate gbias and counts using the extant value of gmult. Similar commands follow for evaluating gmult and gbias, with most optimization combinations possible in the B,G,C command hierarchy (see appendix 2 for a list).

Many of these commands will output the log of the likelihood and/or inquire if a reduced likelihood source test is desired (i.e., compare lnL for an estimate of ck and gb to that of an estimate of gb only with ck=0, both estimates done with gm fixed, or conversely compare LnL for an estimate of ck and gm to an estimate of gm with ck=0, both done with gb fixed). In addition, some of these commands have counterparts that allow maps of each parameter to be produced (the M command set) as described in section 6.4. One command of note is GBC, which evaluates counts, gmult, and gbias simultaneously. The wide range of commands exist to provide as much flexibility as possible. for most EGRET analyses, the L commands pertain.

6.3  Likelihood Testing

Once the model parameter estimation is satisfactory, the likelihood ratio test is performed to determine the significance, or alternatively, the upper limits for the location of interest. Some of the commands noted above perform the likelihood ratio test but the Lxx command set provides for a complete set of tests of both source existence and source position for both a previously undetected source as well as sources from input catalogs.

In terms of source existence, the most elementary functions and a possible starting point, are the LL and L commands which will perform a counts estimate for the extant and input position, respectively. These commands will produce full estimates of all model parameters, calculate the likelihood test statistic (i.e., compare lnL for an estimate of ck, gm, and gb to that for an estimate of gm and gb with ck = 0), calculate the gamma-ray flux, and report other relevant information. These can be followed up with a likelihood position analysis, the most fundamental of which would use the LE command which will find the location of the maximum test statistic as well as the centroid of the 95% confidence source location region, and an elliptical fit to the border of this region and plot it to the selected graphics device. LE makes the position of the maximum TS the source position for unidentified sources. In general, for unidentified sources, it is an open question as to which coordinates should be taken as the source position. Currently, the centroid is the most commonly used value.

Processing of multiple sources can be accomplished through the LM command. LM provides for an automatic iterative solution for multiple PSFs - those in the OTHER PSF map which are marked as active and can be produced or adjusted by the PA command mentioned above. The LM command is very useful when multiple sources exist in the same region. Since the counts assigned to a passive source in the OTHER PSF model affects the analysis of other sources, a simultaneous solution is required. The LM command achieves this through iterative analysis of sources in the OTHER PSF list. The LM command also allows for position analysis (invoking the same analysis as the LE command for each source), or input positions may be assumed. The output of LM to a disk file (with position analysis) is intended to be used in the creation of the EGRET team phase 1 catalog. Also, a command LC exists for an automatic analysis of all sources on the map in an input catalog. It creates a catalog of EGRET fluxes (and upper limits). The flux for sources with TS > TS_max is presented as a detection only. The flux for sources with TS_min < TS < TS_max is presented both as a detection and as an upper limit. The flux for sources with TS < TS_min is presented as an upper limit only (both TS_min and TS_max are adjustable with the PMT command). The default value TS_max is 25. The default value TS_min is 1. The upper limits are determined with a golden search for the Counts value which results in lnL falling by a value of delta_lnL from the value at Counts=0. The default value of delta_lnL is 3.92/2 corresponding to a 95% confidence level. The value of delta_lnL may be changed with the PML command.

Finally, the LP command set are also useful for estimating source positions. For instance, LPS invokes the AMOEBA search for the position which gives the maximum TS for a source. The LPO command will do this amoeba position optimization iteratively on each source marked in the OTHER PSF list until source positions stabilize. This is useful since the position assumed for a source can affect the position indicated for a nearby source. The MF command explained in the next section is also useful for position analysis.

After application of the above commands, the analysis should be fairly complete. Experimentation with the full command list shown in Appendix B is the best way to find the commands of value to each investigator. The next section describes some of the ways in which LIKE outputs the results of the analysis.

6.4  LIKE Output

Many LIKE commands routinely produce output graphics which can be viewed at the console or written to Postscript files by the PGPLOT subroutines (the default graphics device being XWINDOW). This section concentrates on those commands specifically designed to output results to files to save or to be further analyzed with other software. These commands fall in the Mx and Oxx hierarchy. This section describes a few of these options.

The Mx commands create FITS maps of some combination of the likelihood test statistic, g_m, g_b and counts. The maps written by these command are the same size as the ROI. Maps are written to NEWMAPFILE, GBIASFILE, or LMAPFILE depending on which option is used. These file definitions can be adjusted in the CTL file as explained in Appendix 3 or changed with the PMM command. Header information comes from the analysis or from the CTL file (i.e., site which can be changed freely). IF LMAPFILE, , GBIASFILE, or NEWMAPFILE are either '', or 'null', no map will be written. If the specified file exists, you will be asked if you want to write over it.

The MF command allows for the choice between solving for Gmult, Gbias, and Counts at each pixel in fine map, or solving for them once in the fine map center (which is much faster, and adequate unless the fine map is large or expected to have structure other than that of a single PSF at the center of the fine map). The map produced is the likelihood TS in the source region. The EGRET SKYMAP program can then be used to display the position confidence contours for these maps as discussed in section 5.2.3. The MS command results in 3 maps being written: TS in LMAPFILE, Gmult in NEWMAPFILE, and Gbias in GBIASFILE. The MG command results in 2 maps being written: counts in NEWMAPFILE, and Gmult in GBIASFILE. Finally, the MH (Map High Latitude) command has been added which calculates a TS for each pixel with a fixed Gmult. 3 maps are written: TS in LMAPFILE, Counts in NEWMAPFILE, and Gbias in GBIASFILE. Normally the MS command pertains. It results in each pixel being successively examined for the existence of a source centered on the pixel. The evidence (or TS) is saved in the LMAPFILE. Examination of the LMAP after the MS command finishes indicates possible high-energy gamma-ray sources.

The Oxx commands typically write maps or produce plots associated with the defined ROI and unlike the Mx commands include flux as a possible parameter. Exceptions are the OMG, OME, and OMC commands which simply write input maps (although PA can be used to include the ''other PSF'' map to the diffuse model before writing with the OMGx command). Also of note is the OMRx command which writes the residual - the difference between the data and the model. For OMG, the residual is normalized by the estimate of the uncertainty. On the other hand, OMRF writes the flux residual while OMRC writes the counts residual. Because of Poisson fluctuations, the OMR maps usually need to be smoothed with the SKYMAP program (e.g., a gaussian with s = 0.75 degrees) before features are apparent. The OMR function can be used to evaluate the quality of the diffuse model, and to see if correct assumptions about sources in the region have been made.

Finally, the OP profile command provides a tool similar to OMR to obtain a subjective evaluation of the quality of the fit of the model to the data. The data is compared to the model for the ROI. All 3 components of the Model are shown. The active PSF must be first added to the other PSF model in order to be shown in the profile. The data and model are summed in one dimension (nominally the smallest dimension of the ROI) and displayed as a function of the other dimension with counts displayed on the ordinate. The OP profile also calculates reduced chi-squared for the data/model comparison profiles.

7  References

Eadie et al. 1971, Statistical Methods in Experimental Physics.

Bertsch, D. L., Dame, T. M., Fichtel, C. E., Hunter, S. D., Sreekumar, P., Stacy, J. G., & Thaddeus, P., 1993, ApJ, in press.

Cash, W., 1979, ApJ, 228, 939.

Kendall, M.G., & Stuart, A. 1961, The Advanced Theory of Statistics, Volume 2, Haffner, New York.

Mattox, J.R., et al. 1994, in preparation for submission to ApJ.

Pollock, A.M., et al., 1981, AA, 94, 116.

Pollock, A.M., et al., 1985, AA 146, 352.

Thompson, D.J. et al., 1993, Ap. J. Supp., in press.

Willis, T., Mattox, J.R., et al., 1993, Proceedings of the second Compton Observatory Symposium, College Park, MD, Sept 20-22, 1993.

Press et al., 1986 Numerical Recipes.

8  Appendix A. A Sample LIKE Session

A collection of console sessions for exemplary LIKE analyses has been established in the directory /pub/software/egret_like/sessions. These examples are available via anonymous ftp from the node antwrp.gsfc.nasa.gov (128.183.4.16) and should help in understanding the methodology involved in accomplishing most LIKE tasks. As an example, the following example shows an attempt to detect the BL Lac object MRK 421 from the viewing period 5.0 data available through the CGRO archive. In the rest of this appendix, bold text implies comments while standard text is LIKE output.

bronco[23]like
 LIKE  - version 4.6 
 Which ctl file? (return for CTL in current directory)

 Will use CTLFILE: CTL                                               
 READING CTL CARDS
Enter CMAPFILE name:counts.vp0040.g002.fits
 CMAPFILEcounts.vp0040.g002.fits                           
 MAPFILE
 counts.vp0040.g002.fits
         
Reading file:counts.vp0040.g002.fits                                
 
--------------------------------------------------------------------------------
 
THERE ARE            2 IMAGES with energy ranges:
            1          30 MeV -          100 MeV
            2         100 MeV -        99999 MeV
 WHICH DO YOU WANT? 
2
Reading file:exposr.vp0040.g002.fits                                
 
--------------------------------------------------------------------------------
 
MPE software does not support PSF for this energy range;
  PSF being generated for   100.0000     < E <   10000.00    
 Using CALFILEs with suffix:07
 /analysis1/data/difmaps/cfgas.cel.g002b           
 is not installed. Will read unconvolved map.
Reading file:/analysis1/data/difmaps/gas.                           
 
--------------------------------------------------------------------------------
 
NO: 
 /analysis1/data/difmaps/gas.
         
 Error code (N) found in LIKE    ; message is:
MAPRED: MAPFILE DOESN'T EXIST

 Execution continues.
Enter GMAPFILE name:cfgas.cel.g002b
Reading file:cfgas.cel.g002b                                        
 
--------------------------------------------------------------------------------
 
Multiplying the gasmap by exposure map.
The model has already been convolved.
Input radius for analysis (Ranal, cr for  15.00):
 Setting ROI to entire map:
 RA.   130.2500     to   234.7500    ; Dec.  0.2500000     to
79.75000    
 Setting up logarithm table for likelihood.
 Setting OTHER PSF map to null.
Do you want to be notified if  0.60 < Gmult < 1.30
and 0.00 < Gbias < 2.00 (cr for no)

B,C,E,G,L,M,O,P,Q,S,X,? >>
Setup has completed. We chose the map from viewing period 4.0 of phase 1 and chose the E > 100 MeV map. Note that LIKE could not find the diffuse map in the expected places so an alternative map file was input. We used the default value of RANAL. Let's go right to optimizing all source parameters for the Catalog position of MRK 421

B,C,E,G,L,M,O,P,Q,S,X,? >>gbc
Input counts test point:   RA (J2000) and DEC
 (A to abort, cr for  182.250  39.750):166.110  38.210
Input Source Name (cr for                   ):MRK 421
Gmult    1.496402 +/-    0.084828
Gbias    0.738028 +/-    0.051214
Counts      50.54 +/-      12.18
Corresponding flux   14.448708 +/-    3.480921 10^-8 cm^-2 s^-1
Log likelihood:     -2355.10
So the optimization went smoothly and flux is apparent from MRK 421. Let'ss get the maximum likelihood test statistic for the current catalog position

B,C,E,G,L,M,O,P,Q,S,X,? >>ll
 --------------------------------------------------------------
         
 Estimate of Gmult and Gbias with counts=0:
Gmult    1.851476 +/-    0.085468
Gbias     0.60170 +/-     0.05158
Log likelihood:     -2368.35
 --------------------------------------------------------------
         
 Simultaneous estimate of Gmult, Gbias and Counts:
Gmult    1.496393 +/-    0.084829
Gbias    0.738041 +/-    0.051214
Counts      50.54 +/-      12.18
Log likelihood:     -2355.10
Test statistic:    26.51
 --------------------------------------------------------------
         
  Name               RA     DEC    sqrt(TS)  Flux+/- 1 sigma    (U.L.)
     Cnts   +/- 1 sigma  (U.L.)    Gmult    Gbias  Ranal  Asp. EXP
  lnL
MRK 421             166.11  38.21    5.1     14.449    3.481
     50.54      12.18              1.496    0.738  15.0  11.0 34.98
-2355.
A test statistic of 26.51 corresponds to a probability of less than 10-6 of the ck=0 model being true (pre-trials). The square root of the TS is the gaussian s so MRK 421 is detected at the 5.1 s level. Let's see what we can find out about the source position with the LE command.

B,C,E,G,L,M,O,P,Q,S,X,? >>le
Input Finemap center       RA (J2000) and DEC
 (A to abort, cr for  166.110  38.210):
 Do you want to calculate Counts, Gmult, and Gbias once at fine map
center, 
 or at each point - which takes ~10 times longer but
  is more accurate at edges (cr for once at map center)?
When in doubt, take the default!

Enter number rows & columns in fine map (cr for 10):
Enter percentage confidence for source location region
(tabulated values: 50, 68, 95, 99; cr for 95%):
 Obtain fine map centered on   166.1100       38.21000    
 Solving for lnL with Counts=0 for entire fine map.
Gmult    1.851461 +/-    0.085468
Gbias     0.60171 +/-     0.05158
 Solving for Counts, Gmult, and Gbias to be used for entire fine map.
Gmult    1.496397 +/-    0.084829
Gbias    0.738039 +/-    0.051214
Counts      50.54 +/-      12.18
 Try finemap full width=   6.449052    
 Maximum TS =   26.50879    ; found at L=   166.1100    ; B=
38.21000    
Map is too large; fraction of pixels in contour is   0.07
 Obtain fine map centered on   166.5399       38.13834    
 Solving for lnL with Counts=0 for entire fine map.
Gmult    1.901675 +/-    0.085219
Gbias     0.57720 +/-     0.05137
 Solving for Counts, Gmult, and Gbias to be used for entire fine map.
Gmult    1.494991 +/-    0.084595
Gbias    0.745041 +/-    0.051015
Counts      50.77 +/-      12.27
 New finemap center=   166.5399       38.13834    
 Try finemap full width=   3.932868    
 Maximum TS =   27.31836    ; found at L=   166.1466    ; B=
38.53163    
Map is too large; fraction of pixels in contour is   0.12
 Try finemap full width=   3.096331    
 Maximum TS =   27.37305    ; found at L=   166.2303    ; B=
38.44798    
 Mean position within error region   166.3851       38.35767    
 Characteristic radius (  0.7589846    ) degrees yields RMS deviation
   5.5977535E-02
 elipse fit (in degrees): major axis a=  0.8072745    ; minor axis b=
   0.6938947    ; position angle phi=   155.6362    ; yields RMS
deviation
   3.9397180E-02 degrees
 The circle with radius Rb=sqrt(a*b)=  0.7484407    
  degrees yields RMS deviation  5.5825576E-02
Plot fit for elipse or circle 
(cr for elipse, Q to return to main menu)
graphics device? (return for /XWIN     ):
Type <RETURN> to continue: Plot fit for elipse or circle 
(cr for elipse, Q to return to main menu)q

The plot of the Likelihood statistic in the region
of MRK 421 for position analysis

Figure 8.1 The plot of the Likelihood statistic in the region of MRK 421 for position analysis. larger
So after a couple of false starts, LIKE produced the goods and reported the location of maximum TS and the centroid of the 95% contour. Figure 8.1 shows the resulting xwindows map. Looks good so now lets try to get a profile plot of the region. First lets's change the ROI to define a swath of sky in declination along the direction of MRK 421.

B,C,E,G,L,M,O,P,Q,S,X,? >>s
 The Region Of Interest (ROI) is a selected rectangle within the map.
 It is the region profiled (command OP), and the region for which
 each pixel is examined for a point source (with the MS command in
LIKE).
 Note that the ROI does not influence which data is used by LIKE for a
point
 source analysis for a specific point. This is always the region
within the
 map within Ranal of the pixel containing the point.

 Current ROI is:
 RA.   130.2500     to   234.7500    ; Dec.  0.2500000     to
79.75000    

ROI_ADJ:Q; T; center; ra1,ra2,d1,d2; .; or ?>>62 70 0 80
 Error code (O) found in ROIADJ  ; message is:
Requested ROI not entirely within the map. Border(s) adjstd.
OOPS! try again

 Execution continues.
 RA.   130.2500     to   130.2500    ; Dec.  0.2500000     to
79.75000    

ROI_ADJ:Q; T; center; ra1,ra2,d1,d2; .; or ?>>

B,C,E,G,L,M,O,P,Q,S,X,? >>s
 Current ROI is:
 RA.   130.2500     to   130.2500    ; Dec.  0.2500000     to
79.75000    

ROI_ADJ:Q; T; center; ra1,ra2,d1,d2; .; or ?>>162 170 0 80
 Error code (O) found in ROIADJ  ; message is:
Requested ROI not entirely within the map. Border(s) adjstd.

 Execution continues.
 RA.   162.2500     to   170.2500    ; Dec.  0.2500000     to
79.75000    

ROI_ADJ:Q; T; center; ra1,ra2,d1,d2; .; or ?>>

And let's also put MRK 421 in the ``other PSF'' map so that the profile will include the MRK 421 PSF

B,C,E,G,L,M,O,P,Q,S,X,? >>pa
You have elected to alter the other PSF map.
The current other PSF map contains:
 No sources in other PSF map.

PSF_ADJ:A,D,I,Q,?>>a
The current other PSF map contains:
   #   NAME                  POSITION       CNTS   Sp. I.   Active
  1   MRK 421            166.11  38.21      50.8  2.00       1.

PSF_ADJ:A,D,I,M[n],O,Q,R[n],T,X[n],.,?>>q

The profile plot showing the PSF at the position of
MRK 421

Figure 8.2 The profile plot showing the PSF at the position of MRK 421. larger
now use the OP command to plot the profile of the ROI. The results are shown in Figure 8.2.

B,C,E,G,L,M,O,P,Q,S,X,? >>op
Will plot ROI: 162.25 to 170.25; and   0.25 to  79.75
OK? (cr to proceed with plot, P to adjust plot parameters, or R to
return to main menu)
Type <RETURN> to continue: Will plot ROI: 162.25 to 170.25; and   
0.25
to  79.75
OK? (cr to proceed with plot, P to adjust plot parameters, or R to
return to main menu)p
graphics device? (return for /XWIN     ):/ps
Looks like you want to graph over latitude    - OK (cr)?
How many pixels to combine along abscissa (cr for   1) ?
Y axis Max value? (cr -> use data)
Will plot ROI: 162.25 to 170.25; and   0.25 to  79.75
OK? (cr to proceed with plot, P to adjust plot parameters, or R to
return to main menu)

Looks good . . . . must be ready to publish . . . .

B,C,E,G,L,M,O,P,Q,S,X,? >>pmp
 publish flag set to  T

9  Appendix B. Commands

These commands are standalone

Bx commands are used to evaluate Gbias.

Cxx commands are used to evaluate source counts at a specific point.

Gxx commands are used to evaluate Gmult.

Lx is for a likelihood ratio test.

Mx commands are for producing the TS.

Oxx commands are for making plots or writing maps.

Px commands adjust parameters.

10  Appendix C. The CTL file

The CTL file is an input file which sets up many of the LIKE parameters and can control several aspects of LIKE I/O. Each line is a control card with the first word in each line specifies the purpose of the line (with REMARK being a comment). The CTL file should rarely be changed and then only with caution. The following example shows a typical CTL file.


HEADER    Annotated control Card File for LIKE.
REMARK    Read as formatted input - parameters must be in correct 
positions.
REMARK    Two documentation cards follow.   
TITLE1    CTL cards for SOURCE:like V 4.7.
TITLE2    This documentation card is not currently used.
REMARK    The following card has no important consequences.
SITE      <SITE  :G> <SYSTEM:U>            <EXETYP:I> 
REMARK    Use EGRET PSF.
PSM4      
<PSFTYP:EGRETCAL....>                                         
REMARK    EGRET PSF generation parameters. Don't change.
PSM1      <PSFINC :  0.2>                 <PSMSZ  :   81>    
(<=201)
PSM2      <CLAT   :  0.0>                 <PGRID  :    
1>               
SELECT2   <ITRMOD :   74> <THETA  :123000000> 
<PHI:123> 
REMARK    Calibration result file specification.
CALSPC    <Event class (curve index):1> <TASC coincidence 
mode:1>
REMARK    curve index 1 is all energy classes (2=A,3=B,4=C,5=A+C)
REMARK    TASC coincidence mode 1 is TASC in coincidence (0 is out)
REMARK    LONGBINS=1 -> map parameters will be determined from counts 
map
MAPMTX    <BINSIZE: 0.00> <LONGBINS:  1>
SELECT1   <LOWLIM :   00> <HIGHLIM:   00> <GAMMA  : 
2.00> <RVAL :20.0>
REMARK    Counts map specification.
REMARK    'ask' will cause the program to ask the operator for the Counts 
map.
FILENAM   <CMAPFILE  :ask
REMARK    Exposure map specification.  
REMARK    if EMAPFILE = 'as cnts', the map in the directory with the same 
REMARK    number as the Counts map will be read.
FILENAM   <EMAPFILE  :as cnts
REMARK    Background map specification.
REMARK    if GMAPFILE = null, no map will be read and the GMAP will be set 
to 0
REMARK    'ask' will cause the program to ask the operator for the gas map.
REMARK    'standard' will cause the program to read the standard diffuse 
map.
FILENAM   <GMAPFILE  :standard
FILENAM   <RESULTFILE:console                              
REMARK    Output map specification.
FILENAM   <LMAPFILE  :Lmap 
FILENAM   <GBIASFILE :Bmap
FILENAM   <NEWMAPFILE:newmap
ENDMAP                                                                  
ENDALL                                                                  


File translated from TEX by TTH, version 1.56.