Pointlike Tutorial: Part I - Run Pointlike Interactively
Pointlike Tutorial:
Also see: Examples: Advanced Reading: Matthew Kerr's PhD. thesis describes lots of details about how pointlike was developed; chapters 3 and 4 are particulary recommended. |
Implemented using series of python objects, including numpy, scipy, and matplotlib (all of which are distributed with GLAST_EXT), pointlike is a framework for doing maximum likelihood source analysis of LAT data; it is distributed with the ScienceTools.
With pointlike, you can perform a spectral analysis; create an ROI; interactively modify an ROI, save an ROI to a file; create TSMaps; plot spatial maps; and associate LAT sources with multiwavelength counterparts with an interactive interface with source-by-source associations (rather than with full catalogs at once).
When to use pointlike. If you think of performing a likelihood analysis as basically consisting of 3 parts, i.e., 1. getting the data; 2. creating a model; and 3. fitting the model, pointlike is a valuable tool to use when creating the model.
Note: gtlike is considered the "gold standard" and is used for step 3, when performing a final fit prior to publication. As always when selecting among different tools for different tasks, there are tradeoffs. For example, one of pointlike's advantages is its speed while gtlike is slower but achieves slightly greater accuracy.
Recommended: It is strongly recommended that you browse through the following examples before proceeding:
Assumptions
For the purposes of this tutorial, it is assumed that you are:
- running on a SLAC Central Linux RHEL5-64 machine in a bash shell
from within your data directory.
- you have performed the Pointlike: Prerequisite Steps and you:
- are using ScienceTools-09-23-02, or higher and
- have access to the FTOOLS and the GLAST_EXT libraries.
Note: When this was written on 05/23/2011, ScienceTools-09-23-01 was the highest release version and did not incorporate the latest changes to pointlike.
Only with 09-23-02, or higher can you:
* get the roi.tsmap function
* get the roi.plot_sed function
09-23-02 will also eliminate the need
for the CALDB workaround, which is
currently included in the pointlike
example_setup script.
To perform a spectral analysis:
- From the prompt enter: ipython
- Then, from the interactive python prompts, enter:
[1]: from skymaps import SkyDir
[2]: from uw.like.pointspec import DataSpecification, SpectralAnalysis
[3]: ds=DataSpecification(binfile="binned_4bpd.fits", ltcube="expCube.fits", ft1files="pyLike_6mosData_binned-ft1.fits", ft2files="pyLike_6mosData_binned-ft2-30s.fits")
Note: The "binned_4bpd.fits" file is generated automatically, and if you did not generated an exposure cube using gtltcube, "expCube.fits" will also be generated
"on-the-fly".[4]: center=SkyDir(086.1111,-38.1838,SkyDir.GALACTIC)
Note: 3C454.3 is EMS 1524; i.e., the center of the ft1 file can be:
center=SkyDir(343.6566,16.1494,SkyDir.EQUATORIAL); or the
direction of 3C 454.3, in which case:
center=SkyDir(086.1111,-38.1838,SkyDir.GALACTIC)[5]: sa=SpectralAnalysis(ds,binsperdec=4,emin=100,emax=100000, irf="P6_V3_DIFFUSE",roi_dir=center,maxROI=10,minROI=10)
Note: As shown in the readback below, the software did not recognize the LivetimeCube file that was generated by gtltcube (3C454_expCube.fits), but generated another one (3C454-expCube.fits) "on-the-fly". At the time this was written, 5/24/2011, no reason for this was readily available, but it is expected that this will not occure when running ScienceTools-09-23-02, or higher.
LivetimeCube file 3C454_expCube.fits does not exist: will generate it from the ft2 files LivetimeCube: Filling 1484/49152 pixels .....................! 1.03264e+10 loading FT2 file pyLike_6mosData_binned-ft2-30s.fits .....set Data theta cut at 66.4 deg loading 1 FT1 file(s) pyLike_6mosData_binned-ft1.fits...pyLike_6mosData_binned-ft1.fits Reading from FT1 file list. Cuts: zenith theta and instrument theta (min,max): 105, 0, 66.4 pyLike_6mosData_binned-ft1.fits: found interval 243756000-259308000 photons found: 71172, kept: 52951 (total: 52951) .....saving binfile binned_4bpd.fits for subsequent use .....loading binfile binned_4bpd.fits ...Creating Bands, pixels from file binned_4bpd.fits, tables BANDS, PIXELS GTI interval: 243756000-259308000, on time: 1.27058e+07 s found 28 bands, energies 0-1000000 MeV |
[6]: roi=sa.roi_from_xml(roi_dir=center,xmlfile="example_3C454-srcModel.xml")
.....setting up point sources (12 in ROI)... done!
.....setting up diffuse/extended backgrounds for 24 bands...
.......... EG_v02
.......... GAL_v02 |
[7]: print roi
------> print(roi)
Warning: The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze
the integrand in order to determine the difficulties. If the position of a
local difficulty can be determined (singularity, discontinuity) one will
probably gain from splitting up the interval and calling the integrator
on the subranges. Perhaps a special-purpose integrator should be used. |
# Note: When performing a spectral fit, it is always recommended that you use use_gradient; the following command is like running gtlike:
[8]: roi.fit(use_gradient=True)
.....performing likelihood maximization... Did not converge on first gradient iteration. |
[9]: print roi
------> print(roi)
Warning: The integral is probably divergent, or slowly convergent.
Warning: The integral is probably divergent, or slowly convergent.
Warning: The integral is probably divergent, or slowly convergent.
======== POINT SOURCE FITS ============ |
Caution! The following step cannot be performed unless using ScienceTools-09-23-02, or higher.
[10]: roi.plot_sed(which='EMS1524',outdir-'sed_3C454_before.png)
Note that a simple powerlaw spectral model is not a good fit to the data!
Replace with a spectral model which is curved:
[11]: from uw.like.Models import LogParabola
[12]: roi.modify(which='EMS1524',model=LogParabola())
[13]: roi.fit(use_gradient=True)
.....performing likelihood maximization... Function value at minimum: 22268.153
Attempting to invert full hessian...
Found NaN in covariance matrix!
Skipping full Hessian inversion, trying point source parameter subset...
Found NaN in covariance matrix!
Error in calculating and inverting hessian.
Out[13]: -22268.153393227374 |
Caution! The following step cannot be performed unless using ScienceTools-09-23-02, or higher.
[14]: roi.plot_sed(which='EMS1524',outdir='sed_3C454_after.png')
Notice that the log parabola sectral model fits much better! This is like running gtlike for several smaller energy bins:
# Fit the position of the source. Note that it moved a little bit, so
# it wasn't originally in quite the best position! ??????????????????????????
The following is like running gtsrcfind:
[15]: roi.localize(which='EMS1524')
Localizing source EMS1524, tolerance=1.0e-03... moved delta ra dec a b qual 0.0000 0.0000 343.4960 16.1523 0.0062 0.0062 343.5001 16.1474 0.0018 0.0017 17.0029 0.0009 0.0067 343.5009 16.1475 0.0093 0.0089 0.6240 0.0000 0.0068 343.5010 16.1475 0.0095 0.0091 0.0170 TS change: 0.52 Out[15]: (SkyDir(343.501,16.148), 1, 0.0067643667404888213, 0.51611033391964156) |
Make a test statistic map. This is a residual TS maps, which can be used to look for problems in the model. Generally, there are no very large values in the TS map, which means we have a pretty good model of the sky.
The following step is like running gttsmap:
Caution! This step cannot be performed unless using
ScienceTools-09-23-02, or higher.[16]: roi.tsmap(outfile='tsmap.fits', size=5)
Get the TS value of a source. TS is a measure of significance.
[17]: ts=roi.TS(which='EMS1524')
[18]: print 'The TS of EMS1524 is %s' % ts
-------> print('The TS of EMS1524 is %s' % ts)
The TS of EMS1524 is 14239.4032672 |
The following saved file can be run with gtlike to crosscheck the analysis:
[19]: roi.toXML('saved_results.xml')
Warning Found Index=19.6159925396 > 5, maximum allowed value, Setting parameter value to maximum. Warning Found Index=85.5375723678 > 5, maximum allowed value, Setting parameter value to maximum. Warning Found Index=25.0755867167 > 5, maximum allowed value, Setting parameter value to maximum. |
This region file can be opened in ds9 to see your sources:
[20]: roi.toRegion('saved_results.reg')
Create a results file just like gtlike creates:
[21]: roi.toResults('saved_results.dat')
Owned by:Joshua Lande |
|
| Last updated by: Chuck Patterson 05/24/2011 |