Commit 892c9f419f6168faa5c8bd92ba39d0624acccbda
1 parent
45ba4182
Add stacked cube analysis script.
The purpose of the make_cube_analysis.py script is to provide an end-to-end test fur stacked cube analysis. It uses 4 CTA observations to illustrate the usage of stacked cube analysis and to demonstrate that the fitting results are okay.
Showing
1 changed file
with
154 additions
and
0 deletions
examples/make_cube_analysis.py
0 → 100755
| 1 | +#! /usr/bin/env python | ||
| 2 | +# ========================================================================== | ||
| 3 | +# This script illustrates how to perform a stacked cube analysis. | ||
| 4 | +# | ||
| 5 | +# Copyright (C) 2014 Juergen Knoedlseder | ||
| 6 | +# | ||
| 7 | +# This program is free software: you can redistribute it and/or modify | ||
| 8 | +# it under the terms of the GNU General Public License as published by | ||
| 9 | +# the Free Software Foundation, either version 3 of the License, or | ||
| 10 | +# (at your option) any later version. | ||
| 11 | +# | ||
| 12 | +# This program is distributed in the hope that it will be useful, | ||
| 13 | +# but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 14 | +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| 15 | +# GNU General Public License for more details. | ||
| 16 | +# | ||
| 17 | +# You should have received a copy of the GNU General Public License | ||
| 18 | +# along with this program. If not, see <http://www.gnu.org/licenses/>. | ||
| 19 | +# | ||
| 20 | +# ========================================================================== | ||
| 21 | +import ctools | ||
| 22 | +import gammalib | ||
| 23 | + | ||
| 24 | + | ||
| 25 | +# ================= # | ||
| 26 | +# Analysis pipeline # | ||
| 27 | +# ================= # | ||
| 28 | +def pipeline(caldb = "dummy", \ | ||
| 29 | + irf = "cta_dummy_irf", \ | ||
| 30 | + model = "$CTOOLS/share/models/crab.xml"): | ||
| 31 | + """ | ||
| 32 | + Analysis pipeline for a stacked cube analysis. The pipeline writes | ||
| 33 | + intermediate results to FITS files on disk. Also log files of all | ||
| 34 | + individual ctools will be written. | ||
| 35 | + | ||
| 36 | + Keywords: | ||
| 37 | + caldb - Calibration database repository | ||
| 38 | + irf - Name of instrument response | ||
| 39 | + model - Source model | ||
| 40 | + """ | ||
| 41 | + # Set pointing pattern (4 pointings) | ||
| 42 | + pointings = ctools.obsutils.set_obs_patterns("four", offset=2.0) | ||
| 43 | + | ||
| 44 | + # Set CTA observations (default: 1800 sec per pointing) | ||
| 45 | + obs = ctools.obsutils.set_obs_list(pointings, caldb=caldb, irf=irf) | ||
| 46 | + | ||
| 47 | + # Load models | ||
| 48 | + models = gammalib.GModels(model) | ||
| 49 | + obs.models(models) | ||
| 50 | + | ||
| 51 | + # Run ctobssim (simulates 4 event lists) | ||
| 52 | + sim = ctools.ctobssim(obs) | ||
| 53 | + sim.logFileOpen() | ||
| 54 | + sim.run() | ||
| 55 | + | ||
| 56 | + # Run ctbin (builds stacked event cube) | ||
| 57 | + bin = ctools.ctbin(sim.obs()) | ||
| 58 | + bin.logFileOpen() | ||
| 59 | + bin["outfile"].filename("cntmap.fits") | ||
| 60 | + bin["ebinalg"].string("LOG") | ||
| 61 | + bin["emin"].real(0.1) | ||
| 62 | + bin["emax"].real(100.0) | ||
| 63 | + bin["enumbins"].integer(20) | ||
| 64 | + bin["usepnt"].boolean(False) | ||
| 65 | + bin["nxpix"].integer(300) | ||
| 66 | + bin["nypix"].integer(300) | ||
| 67 | + bin["binsz"].real(0.02) | ||
| 68 | + bin["coordsys"].string("CEL") | ||
| 69 | + bin["proj"].string("CAR") | ||
| 70 | + bin["xref"].real(83.6331) | ||
| 71 | + bin["yref"].real(22.0145) | ||
| 72 | + bin.execute() | ||
| 73 | + | ||
| 74 | + # Run ctexpcube (builds exposure cube) | ||
| 75 | + expcube = ctools.ctexpcube(sim.obs()) | ||
| 76 | + expcube.logFileOpen() | ||
| 77 | + expcube["cntmap"].filename("cntmap.fits") | ||
| 78 | + expcube["outfile"].filename("expcube.fits") | ||
| 79 | + expcube["caldb"].string(caldb) | ||
| 80 | + expcube["irf"].string(irf) | ||
| 81 | + expcube.execute() | ||
| 82 | + | ||
| 83 | + # Run ctpsfcube (builds PSF cube) | ||
| 84 | + psfcube = ctools.ctpsfcube(sim.obs()) | ||
| 85 | + psfcube.logFileOpen() | ||
| 86 | + psfcube["cntmap"].filename("NONE") | ||
| 87 | + psfcube["outfile"].filename("psfcube.fits") | ||
| 88 | + psfcube["caldb"].string(caldb) | ||
| 89 | + psfcube["irf"].string(irf) | ||
| 90 | + psfcube["ebinalg"].string("LOG") | ||
| 91 | + psfcube["emin"].real(0.1) | ||
| 92 | + psfcube["emax"].real(100.0) | ||
| 93 | + psfcube["enumbins"].integer(20) | ||
| 94 | + psfcube["nxpix"].integer(5) | ||
| 95 | + psfcube["nypix"].integer(5) | ||
| 96 | + psfcube["binsz"].real(2.0) | ||
| 97 | + psfcube["coordsys"].string("CEL") | ||
| 98 | + psfcube["proj"].string("CAR") | ||
| 99 | + psfcube["xref"].real(83.6331) | ||
| 100 | + psfcube["yref"].real(22.0145) | ||
| 101 | + psfcube["amax"].real(0.3) | ||
| 102 | + psfcube["anumbins"].integer(60) | ||
| 103 | + psfcube.execute() | ||
| 104 | + | ||
| 105 | + # Run ctbkgcube (builds background cube) | ||
| 106 | + bkgcube = ctools.ctbkgcube(sim.obs()) | ||
| 107 | + bkgcube.logFileOpen() | ||
| 108 | + bkgcube["cntmap"].filename("cntmap.fits") | ||
| 109 | + bkgcube["bkgmdl"].filename(model) | ||
| 110 | + bkgcube["outfile"].filename("bkgcube.fits") | ||
| 111 | + bkgcube.execute() | ||
| 112 | + | ||
| 113 | + # Create observation container for stacked cube analysis | ||
| 114 | + obs.clear() | ||
| 115 | + cta = gammalib.GCTAObservation("cntmap.fits", "expcube.fits", "psfcube.fits") | ||
| 116 | + obs.append(cta) | ||
| 117 | + | ||
| 118 | + # Create model for stacked cube analysis | ||
| 119 | + models = gammalib.GModels(model) | ||
| 120 | + models.remove("Background") | ||
| 121 | + spatial = gammalib.GModelSpatialDiffuseCube("bkgcube.fits") | ||
| 122 | + spectral = gammalib.GModelSpectralConst() | ||
| 123 | + background = gammalib.GCTAModelCubeBackground(spatial, spectral) | ||
| 124 | + models.append(background) | ||
| 125 | + models.save("cube_model.xml") | ||
| 126 | + | ||
| 127 | + # Append model to observation container | ||
| 128 | + obs.models(models) | ||
| 129 | + | ||
| 130 | + # Run ctlike | ||
| 131 | + like = ctools.ctlike(obs) | ||
| 132 | + like.logFileOpen() | ||
| 133 | + like["debug"].boolean(True) | ||
| 134 | + like["outmdl"].filename("cube_results.xml") | ||
| 135 | + like.execute() | ||
| 136 | + | ||
| 137 | + # Return | ||
| 138 | + return | ||
| 139 | + | ||
| 140 | + | ||
| 141 | +#==========================# | ||
| 142 | +# Main routine entry point # | ||
| 143 | +#==========================# | ||
| 144 | +if __name__ == '__main__': | ||
| 145 | + """ | ||
| 146 | + """ | ||
| 147 | + # Dump header | ||
| 148 | + print("*************************************") | ||
| 149 | + print("* CTA stacked cube analysis *") | ||
| 150 | + print("*************************************") | ||
| 151 | + | ||
| 152 | + # Run pipeline | ||
| 153 | + pipeline() | ||
| 154 | + | ||
| 0 | \ No newline at end of file | 155 | \ No newline at end of file |