Commit 6418ebd17599adcb809da65422ad1d9ca35e6ac9
Merge branch '1264-obsutils-under-ctools' into integration
Conflicts: scripts/cspull.py scripts/obsutils.py
Showing
12 changed files
with
540 additions
and
534 deletions
.gitignore
| @@ -44,8 +44,8 @@ doc/man | @@ -44,8 +44,8 @@ doc/man | ||
| 44 | src/ctools-setup | 44 | src/ctools-setup |
| 45 | pyext/build/ | 45 | pyext/build/ |
| 46 | pyext/setup.py | 46 | pyext/setup.py |
| 47 | -pyext/ctools.py | ||
| 48 | -pyext/ctools_wrap.cpp | 47 | +pyext/ctools/ctools.py |
| 48 | +pyext/ctools/ctools_wrap.cpp | ||
| 49 | src/ctbin/ctbin | 49 | src/ctbin/ctbin |
| 50 | src/ctlike/ctlike | 50 | src/ctlike/ctlike |
| 51 | src/ctmodel/ctmodel | 51 | src/ctmodel/ctmodel |
configure.ac
| @@ -630,7 +630,7 @@ else | @@ -630,7 +630,7 @@ else | ||
| 630 | 630 | ||
| 631 | # Do we have wrappers (ctools_wrap.cpp and ctools.py)? | 631 | # Do we have wrappers (ctools_wrap.cpp and ctools.py)? |
| 632 | has_wrappers="yes" | 632 | has_wrappers="yes" |
| 633 | - AC_CHECK_FILES([pyext/ctools_wrap.cpp pyext/ctools.py],, | 633 | + AC_CHECK_FILES([pyext/ctools/ctools_wrap.cpp pyext/ctools/ctools.py],, |
| 634 | [has_wrappers="no"]) | 634 | [has_wrappers="no"]) |
| 635 | if test "x$has_wrappers" = "xno"; then | 635 | if test "x$has_wrappers" = "xno"; then |
| 636 | AC_MSG_WARN([Python wrapper missing. Requires swig for wrapper generation.]) | 636 | AC_MSG_WARN([Python wrapper missing. Requires swig for wrapper generation.]) |
| @@ -669,7 +669,7 @@ if test "x$enable_python" = "xyes"; then | @@ -669,7 +669,7 @@ if test "x$enable_python" = "xyes"; then | ||
| 669 | 669 | ||
| 670 | else | 670 | else |
| 671 | if test "x$has_python_wrap" = "xno" -a "x$has_swig" = "xno"; then | 671 | if test "x$has_python_wrap" = "xno" -a "x$has_swig" = "xno"; then |
| 672 | - AC_MSG_WARN([Neither pyext/ctools_wrap.cpp and pyext/ctools.py nor swig has not been found]) | 672 | + AC_MSG_WARN([Neither pyext/ctools/ctools_wrap.cpp and pyext/ctools/ctools.py nor swig has not been found]) |
| 673 | AC_MSG_WARN([Python binding for ctools cannot be built]) | 673 | AC_MSG_WARN([Python binding for ctools cannot be built]) |
| 674 | AC_MSG_WARN([Either use a ctools distribution with ctools_wrap.cpp and ctools.py shipped or install swig on your system]) | 674 | AC_MSG_WARN([Either use a ctools distribution with ctools_wrap.cpp and ctools.py shipped or install swig on your system]) |
| 675 | AC_MSG_WARN([swig can be obtained from http://www.swig.org/]) | 675 | AC_MSG_WARN([swig can be obtained from http://www.swig.org/]) |
pyext/Makefile.am
| @@ -30,7 +30,7 @@ | @@ -30,7 +30,7 @@ | ||
| 30 | if WITH_PYTHON | 30 | if WITH_PYTHON |
| 31 | 31 | ||
| 32 | # Set SWIG interface filename | 32 | # Set SWIG interface filename |
| 33 | -SWIG_INTERFACE = ctools.i | 33 | +SWIG_INTERFACE = ctools/ctools.i |
| 34 | 34 | ||
| 35 | # Dependencies. Obtain the list of dependencies by executing ls. | 35 | # Dependencies. Obtain the list of dependencies by executing ls. |
| 36 | DEPFILES = $(shell ls *.Pi 2> /dev/null) | 36 | DEPFILES = $(shell ls *.Pi 2> /dev/null) |
| @@ -39,7 +39,7 @@ include $(DEPFILES) | @@ -39,7 +39,7 @@ include $(DEPFILES) | ||
| 39 | # Define the targets that should be built before anything else. | 39 | # Define the targets that should be built before anything else. |
| 40 | # These are honored by 'make all', 'make check' and 'make install'. | 40 | # These are honored by 'make all', 'make check' and 'make install'. |
| 41 | if HAVE_SWIG | 41 | if HAVE_SWIG |
| 42 | - BUILT_SOURCES = ctools_wrap.cpp ctools.py | 42 | + BUILT_SOURCES = ctools/ctools_wrap.cpp ctools/ctools.py |
| 43 | endif | 43 | endif |
| 44 | 44 | ||
| 45 | # Set files to be removed in clean. We remove only the SWIG interface files | 45 | # Set files to be removed in clean. We remove only the SWIG interface files |
| @@ -59,7 +59,7 @@ all-local: build | @@ -59,7 +59,7 @@ all-local: build | ||
| 59 | # of this specific directory. So we build one instance for checking, the other | 59 | # of this specific directory. So we build one instance for checking, the other |
| 60 | # for installation. If we could get the name built by distutils we could avoid | 60 | # for installation. If we could get the name built by distutils we could avoid |
| 61 | # this. | 61 | # this. |
| 62 | -build: ctools_wrap.cpp ctools.py | 62 | +build: ctools/ctools_wrap.cpp ctools/ctools.py |
| 63 | $(PYTHON) setup.py build | 63 | $(PYTHON) setup.py build |
| 64 | 64 | ||
| 65 | install-exec-local: | 65 | install-exec-local: |
| @@ -70,9 +70,9 @@ clean-local: | @@ -70,9 +70,9 @@ clean-local: | ||
| 70 | 70 | ||
| 71 | # Rule for building Python wrappers using SWIG. | 71 | # Rule for building Python wrappers using SWIG. |
| 72 | if HAVE_SWIG | 72 | if HAVE_SWIG |
| 73 | -ctools.py: ctools_wrap.cpp | ||
| 74 | -ctools_wrap.cpp: $(SWIG_INTERFACE) | ||
| 75 | - if $(SWIGCOMPILE) -MMD -MF "ctools.Tpi" -o ctools_wrap.cpp $<; \ | 73 | +ctools/ctools.py: ctools/ctools_wrap.cpp |
| 74 | +ctools/ctools_wrap.cpp: $(SWIG_INTERFACE) | ||
| 75 | + if $(SWIGCOMPILE) -MMD -MF "ctools.Tpi" -o ctools/ctools_wrap.cpp $<; \ | ||
| 76 | then mv -f "ctools.Tpi" "ctools.Pi"; else rm -f "ctools.Tpi"; exit 1; fi | 76 | then mv -f "ctools.Tpi" "ctools.Pi"; else rm -f "ctools.Tpi"; exit 1; fi |
| 77 | endif # HAVE_SWIG | 77 | endif # HAVE_SWIG |
| 78 | 78 |
pyext/ctools/__init__.py
0 → 100644
pyext/ctools.i renamed to pyext/ctools/ctools.i
pyext/ctools/obsutils.py
0 → 100755
| 1 | +# ========================================================================== | ||
| 2 | +# This script provides a number of functions that are useful for handling | ||
| 3 | +# CTA observations. | ||
| 4 | +# | ||
| 5 | +# Copyright (C) 2011-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 | +# Simulate observations # | ||
| 27 | +# ===================== # | ||
| 28 | +def sim(obs, log=False, debug=False, edisp=False, seed=0, nbins=0, | ||
| 29 | + binsz=0.05, npix=200, proj="TAN", coord="GAL"): | ||
| 30 | + """ | ||
| 31 | + Simulate events for all observations in the container. | ||
| 32 | + | ||
| 33 | + Parameters: | ||
| 34 | + obs - Observation container | ||
| 35 | + Keywords: | ||
| 36 | + log - Create log file(s) | ||
| 37 | + debug - Create console dump? | ||
| 38 | + edisp - Apply energy dispersion? | ||
| 39 | + seed - Seed value for simulations (default: 0) | ||
| 40 | + nbins - Number of energy bins (default: 0=unbinned) | ||
| 41 | + binsz - Pixel size for binned simulation (deg/pixel) | ||
| 42 | + npix - Number of pixels in X and Y for binned simulation | ||
| 43 | + """ | ||
| 44 | + # Allocate ctobssim application and set parameters | ||
| 45 | + sim = ctools.ctobssim(obs) | ||
| 46 | + sim["seed"].integer(seed) | ||
| 47 | + | ||
| 48 | + # Optionally open the log file | ||
| 49 | + if log: | ||
| 50 | + sim.logFileOpen() | ||
| 51 | + | ||
| 52 | + # Optionally switch-on debugging model | ||
| 53 | + if debug: | ||
| 54 | + sim["debug"].boolean(True) | ||
| 55 | + | ||
| 56 | + # Run ctobssim application. This will loop over all observations in the | ||
| 57 | + # container and simulation the events for each observation. Note that | ||
| 58 | + # events are not added together, they still apply to each observation | ||
| 59 | + # separately. | ||
| 60 | + sim.run() | ||
| 61 | + | ||
| 62 | + # Binned option? | ||
| 63 | + if nbins > 0: | ||
| 64 | + | ||
| 65 | + # Determine common energy boundaries for observations | ||
| 66 | + emin = None | ||
| 67 | + emax = None | ||
| 68 | + for run in sim.obs(): | ||
| 69 | + run_emin = run.events().ebounds().emin().TeV() | ||
| 70 | + run_emax = run.events().ebounds().emax().TeV() | ||
| 71 | + if emin == None: | ||
| 72 | + emin = run_emin | ||
| 73 | + elif run_emin > emin: | ||
| 74 | + emin = run_emin | ||
| 75 | + if emax == None: | ||
| 76 | + emax = run_emax | ||
| 77 | + elif run_emax > emax: | ||
| 78 | + emax = run_emax | ||
| 79 | + | ||
| 80 | + # Allocate ctbin application and set parameters | ||
| 81 | + bin = ctools.ctbin(sim.obs()) | ||
| 82 | + bin["ebinalg"].string("LOG") | ||
| 83 | + bin["emin"].real(emin) | ||
| 84 | + bin["emax"].real(emax) | ||
| 85 | + bin["enumbins"].integer(nbins) | ||
| 86 | + bin["usepnt"].boolean(True) # Use pointing for map centre | ||
| 87 | + bin["nxpix"].integer(npix) | ||
| 88 | + bin["nypix"].integer(npix) | ||
| 89 | + bin["binsz"].real(binsz) | ||
| 90 | + bin["coordsys"].string(coord) | ||
| 91 | + bin["proj"].string(proj) | ||
| 92 | + | ||
| 93 | + # Optionally open the log file | ||
| 94 | + if log: | ||
| 95 | + bin.logFileOpen() | ||
| 96 | + | ||
| 97 | + # Optionally switch-on debugging model | ||
| 98 | + if debug: | ||
| 99 | + bin["debug"].boolean(True) | ||
| 100 | + | ||
| 101 | + # Run ctbin application. This will loop over all observations in | ||
| 102 | + # the container and bin the events in counts maps | ||
| 103 | + bin.run() | ||
| 104 | + | ||
| 105 | + # Make a deep copy of the observation that will be returned | ||
| 106 | + # (the ctbin object will go out of scope one the function is | ||
| 107 | + # left) | ||
| 108 | + obs = bin.obs().copy() | ||
| 109 | + | ||
| 110 | + else: | ||
| 111 | + | ||
| 112 | + # Make a deep copy of the observation that will be returned | ||
| 113 | + # (the ctobssim object will go out of scope one the function is | ||
| 114 | + # left) | ||
| 115 | + obs = sim.obs().copy() | ||
| 116 | + | ||
| 117 | + # Delete the simulation | ||
| 118 | + del sim | ||
| 119 | + | ||
| 120 | + # Return observation container | ||
| 121 | + return obs | ||
| 122 | + | ||
| 123 | + | ||
| 124 | +# ================ # | ||
| 125 | +# Fit observations # | ||
| 126 | +# ================ # | ||
| 127 | +def fit(obs, log=False, debug=False, edisp=False): | ||
| 128 | + """ | ||
| 129 | + Perform maximum likelihood fitting of observations in the container. | ||
| 130 | + | ||
| 131 | + Parameters: | ||
| 132 | + obs - Observation container | ||
| 133 | + Keywords: | ||
| 134 | + log - Create log file(s) | ||
| 135 | + debug - Create screen dump | ||
| 136 | + edisp - Apply energy dispersion? | ||
| 137 | + """ | ||
| 138 | + # Allocate ctlike application | ||
| 139 | + like = ctools.ctlike(obs) | ||
| 140 | + | ||
| 141 | + # Optionally open the log file | ||
| 142 | + if log: | ||
| 143 | + like.logFileOpen() | ||
| 144 | + | ||
| 145 | + # Optionally switch-on debugging model | ||
| 146 | + if debug: | ||
| 147 | + like["debug"].boolean(True) | ||
| 148 | + | ||
| 149 | + # Optionally apply energy dispersion | ||
| 150 | + like["edisp"].boolean(edisp) | ||
| 151 | + | ||
| 152 | + # Run ctlike application. | ||
| 153 | + like.run() | ||
| 154 | + | ||
| 155 | + # Return observations | ||
| 156 | + return like | ||
| 157 | + | ||
| 158 | + | ||
| 159 | +# ================= # | ||
| 160 | +# Create counts map # | ||
| 161 | +# ================= # | ||
| 162 | +def cntmap(obs, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ | ||
| 163 | + binsz=0.05, nxpix=200, nypix=200, \ | ||
| 164 | + outname="cntmap.fits"): | ||
| 165 | + """ | ||
| 166 | + Creates a counts map by combining the events of all observations. | ||
| 167 | + The counts map will be a summed map over all energies. | ||
| 168 | + | ||
| 169 | + Parameters: | ||
| 170 | + obs - Observation container | ||
| 171 | + Keywords: | ||
| 172 | + proj - Projection type (e.g. TAN, CAR, STG, ...) (default: TAN) | ||
| 173 | + coord - Coordinate type (GAL, CEL) (default: GAL) | ||
| 174 | + xval - Reference longitude value [deg] (default: 0.0) | ||
| 175 | + yval - Reference latitude value [deg] (default: 0.0) | ||
| 176 | + binsz - Pixel size [deg/pixel] (default: 0.05) | ||
| 177 | + nxpix - Number of pixels in X direction (default: 200) | ||
| 178 | + nypix - Number of pixels in Y direction (default: 200) | ||
| 179 | + outname - Counts map FITS filename (default: cntmap.fits) | ||
| 180 | + """ | ||
| 181 | + # Allocate counts map | ||
| 182 | + map = gammalib.GSkymap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1) | ||
| 183 | + | ||
| 184 | + # Set maximum pixel number | ||
| 185 | + maxpixel = nxpix * nypix | ||
| 186 | + | ||
| 187 | + # Fill all observations | ||
| 188 | + for run in obs: | ||
| 189 | + | ||
| 190 | + # Loop over all events | ||
| 191 | + for event in run.events(): | ||
| 192 | + | ||
| 193 | + # Determine sky pixel | ||
| 194 | + skydir = gammalib.GCTAInstDir(event.dir()).dir() | ||
| 195 | + pixel = map.dir2inx(skydir) | ||
| 196 | + | ||
| 197 | + # Set pixel | ||
| 198 | + if pixel < maxpixel: | ||
| 199 | + map[pixel] += event.counts() | ||
| 200 | + | ||
| 201 | + # Save sky map. The clobber flag is set to True, so any existing FITS | ||
| 202 | + # file will be overwritten. | ||
| 203 | + map.save(outname, True) | ||
| 204 | + | ||
| 205 | + # Return counts map | ||
| 206 | + return map | ||
| 207 | + | ||
| 208 | + | ||
| 209 | +# ================ # | ||
| 210 | +# Create model map # | ||
| 211 | +# ================ # | ||
| 212 | +def modmap(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ | ||
| 213 | + binsz=0.05, nxpix=200, nypix=200, \ | ||
| 214 | + outname="modmap.fits"): | ||
| 215 | + """ | ||
| 216 | + Make model map for a given reference energy by combining all observations. | ||
| 217 | + The model map will be evaluated for a given reference energy 'eref' and will | ||
| 218 | + be given in units of [counts/(sr MeV s)]. | ||
| 219 | + | ||
| 220 | + Parameters: | ||
| 221 | + obs - Observation container | ||
| 222 | + Keywords: | ||
| 223 | + eref - Reference energy for which model is created [TeV] (default: 0.1) | ||
| 224 | + proj - Projection type (e.g. TAN, CAR, STG, ...) (default: TAN) | ||
| 225 | + coord - Coordinate type (GAL, CEL) (default: GAL) | ||
| 226 | + xval - Reference longitude value [deg] (default: 0.0) | ||
| 227 | + yval - Reference latitude value [deg] (default: 0.0) | ||
| 228 | + binsz - Pixel size [deg/pixel] (default: 0.05) | ||
| 229 | + nxpix - Number of pixels in X direction (default: 200) | ||
| 230 | + nypix - Number of pixels in Y direction (default: 200) | ||
| 231 | + outname - Model map FITS filename (default: modmap.fits) | ||
| 232 | + """ | ||
| 233 | + # Allocate model map | ||
| 234 | + map = gammalib.GSkymap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1) | ||
| 235 | + | ||
| 236 | + # Set reference energy, time and direction. The time is not initialised and is | ||
| 237 | + # in fact not used (as the IRF is assumed to be time independent for now). | ||
| 238 | + # The sky direction is set later using the pixel values. | ||
| 239 | + energy = gammalib.GEnergy() | ||
| 240 | + time = gammalib.GTime() | ||
| 241 | + instdir = gammalib.GCTAInstDir() | ||
| 242 | + energy.TeV(eref) | ||
| 243 | + | ||
| 244 | + # Loop over all map pixels | ||
| 245 | + for pixel in range(map.npix()): | ||
| 246 | + | ||
| 247 | + # Get sky direction | ||
| 248 | + skydir = map.inx2dir(pixel) | ||
| 249 | + instdir.dir(skydir) | ||
| 250 | + | ||
| 251 | + # Create event atom for map pixel | ||
| 252 | + atom = gammalib.GCTAEventAtom() | ||
| 253 | + atom.dir(instdir) | ||
| 254 | + atom.energy(energy) | ||
| 255 | + atom.time(time) | ||
| 256 | + | ||
| 257 | + # Initialise model value | ||
| 258 | + value = 0.0 | ||
| 259 | + | ||
| 260 | + # Loop over all observations | ||
| 261 | + for run in obs: | ||
| 262 | + value += obs.models().eval(atom, run) | ||
| 263 | + | ||
| 264 | + # Set map value | ||
| 265 | + map[pixel] = value | ||
| 266 | + | ||
| 267 | + # Save sky map | ||
| 268 | + map.save(outname, True) | ||
| 269 | + | ||
| 270 | + # Return model map | ||
| 271 | + return map | ||
| 272 | + | ||
| 273 | + | ||
| 274 | +# ===================================== # | ||
| 275 | +# Set one CTA observation (old version) # | ||
| 276 | +# ===================================== # | ||
| 277 | +def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ | ||
| 278 | + emin=0.1, emax=100.0, rad=5.0, \ | ||
| 279 | + irf="cta_dummy_irf", caldb="dummy"): | ||
| 280 | + """ | ||
| 281 | + Obsolete function, use set_obs instead. | ||
| 282 | + """ | ||
| 283 | + # Print warning | ||
| 284 | + print("Warning: obsutils.set is obsolete, use obsutils.set_obs instead.") | ||
| 285 | + | ||
| 286 | + # Call new function | ||
| 287 | + obs_cta = set_obs(pntdir, tstart=tstart, duration=duration, deadc=deadc, \ | ||
| 288 | + emin=emin, emax=emax, rad=rad, \ | ||
| 289 | + irf=irf, caldb=caldb) | ||
| 290 | + | ||
| 291 | + # Return CTA observation | ||
| 292 | + return obs_cta | ||
| 293 | + | ||
| 294 | + | ||
| 295 | +# ======================= # | ||
| 296 | +# Set one CTA observation # | ||
| 297 | +# ======================= # | ||
| 298 | +def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ | ||
| 299 | + emin=0.1, emax=100.0, rad=5.0, \ | ||
| 300 | + irf="cta_dummy_irf", caldb="dummy", id="000000"): | ||
| 301 | + """ | ||
| 302 | + Returns a single CTA observation containing an empty CTA event list. | ||
| 303 | + By looping over this function you can add CTA observations to the | ||
| 304 | + observation container. | ||
| 305 | + | ||
| 306 | + Parameters: | ||
| 307 | + pntdir - Pointing direction [GSkyDir] | ||
| 308 | + Keywords: | ||
| 309 | + tstart - Start time [seconds] (default: 0.0) | ||
| 310 | + duration - Duration of observation [seconds] (default: 1800.0) | ||
| 311 | + deadc - Deadtime correction factor (default: 0.95) | ||
| 312 | + emin - Minimum event energy [TeV] (default: 0.1) | ||
| 313 | + emax - Maximum event energy [TeV] (default: 100.0) | ||
| 314 | + rad - ROI radius used for analysis [deg] (default: 5.0) | ||
| 315 | + irf - Instrument response function (default: cta_dummy_irf) | ||
| 316 | + caldb - Calibration database path (default: "dummy") | ||
| 317 | + """ | ||
| 318 | + # Allocate CTA observation | ||
| 319 | + obs_cta = gammalib.GCTAObservation() | ||
| 320 | + | ||
| 321 | + # Set calibration database | ||
| 322 | + db = gammalib.GCaldb() | ||
| 323 | + if (gammalib.dir_exists(caldb)): | ||
| 324 | + db.rootdir(caldb) | ||
| 325 | + else: | ||
| 326 | + db.open("cta", caldb) | ||
| 327 | + | ||
| 328 | + # Set pointing direction | ||
| 329 | + pnt = gammalib.GCTAPointing() | ||
| 330 | + pnt.dir(pntdir) | ||
| 331 | + obs_cta.pointing(pnt) | ||
| 332 | + | ||
| 333 | + # Set ROI | ||
| 334 | + roi = gammalib.GCTARoi() | ||
| 335 | + instdir = gammalib.GCTAInstDir() | ||
| 336 | + instdir.dir(pntdir) | ||
| 337 | + roi.centre(instdir) | ||
| 338 | + roi.radius(rad) | ||
| 339 | + | ||
| 340 | + # Set GTI | ||
| 341 | + gti = gammalib.GGti() | ||
| 342 | + gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration)) | ||
| 343 | + | ||
| 344 | + # Set energy boundaries | ||
| 345 | + ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \ | ||
| 346 | + gammalib.GEnergy(emax, "TeV")) | ||
| 347 | + | ||
| 348 | + # Allocate event list | ||
| 349 | + events = gammalib.GCTAEventList() | ||
| 350 | + events.roi(roi) | ||
| 351 | + events.gti(gti) | ||
| 352 | + events.ebounds(ebounds) | ||
| 353 | + obs_cta.events(events) | ||
| 354 | + | ||
| 355 | + # Set instrument response | ||
| 356 | + obs_cta.response(irf, db) | ||
| 357 | + | ||
| 358 | + # Set ontime, livetime, and deadtime correction factor | ||
| 359 | + obs_cta.ontime(duration) | ||
| 360 | + obs_cta.livetime(duration*deadc) | ||
| 361 | + obs_cta.deadc(deadc) | ||
| 362 | + obs_cta.id(id) | ||
| 363 | + | ||
| 364 | + # Return CTA observation | ||
| 365 | + return obs_cta | ||
| 366 | + | ||
| 367 | + | ||
| 368 | +# ============================ # | ||
| 369 | +# Set list of CTA observations # | ||
| 370 | +# ============================ # | ||
| 371 | +def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \ | ||
| 372 | + emin=0.1, emax=100.0, rad=5.0, \ | ||
| 373 | + irf="cta_dummy_irf", caldb="dummy"): | ||
| 374 | + """ | ||
| 375 | + Returns an observation container filled with a list of CTA observations. | ||
| 376 | + The list is defined by the obsdeflist parameter which is a dictionnary | ||
| 377 | + containing the mandatory keywords 'ra' and 'dec' that specify the | ||
| 378 | + pointing direction for a given observation. Optional keyword give control | ||
| 379 | + over other observation proposerties, such as duration, deadtime correction, | ||
| 380 | + energy range, etc. If an optional keyword is not specified, the function | ||
| 381 | + keyword is used instead. | ||
| 382 | + | ||
| 383 | + Parameters: | ||
| 384 | + obsdeflist - Observation definition list [{'ra': x.xx, 'dec': x.xx}] | ||
| 385 | + The directory can take the following optional keywords: | ||
| 386 | + - duration: Duration of CTA observation [seconds] | ||
| 387 | + - deadc: Deadtime correction factor | ||
| 388 | + - emin: Minimum event energy [TeV] | ||
| 389 | + - emax: Maximum event energy [TeV] | ||
| 390 | + - rad: ROI radius used for analysis [deg] | ||
| 391 | + - irf: Instrument response function | ||
| 392 | + - caldb: Calibration database path | ||
| 393 | + Optional keywords overwrite keywords specified in the | ||
| 394 | + function call. | ||
| 395 | + Keywords: | ||
| 396 | + tstart - Start time [seconds] (default: 0.0) | ||
| 397 | + duration - Duration of one CTA observation [seconds] (default: 1800.0) | ||
| 398 | + deadc - Deadtime correction factor (default: 0.95) | ||
| 399 | + emin - Minimum event energy [TeV] (default: 0.1) | ||
| 400 | + emax - Maximum event energy [TeV] (default: 100.0) | ||
| 401 | + rad - ROI radius used for analysis [deg] (default: 5.0) | ||
| 402 | + irf - Instrument response function (default: cta_dummy_irf) | ||
| 403 | + caldb - Calibration database path (default: "dummy") | ||
| 404 | + """ | ||
| 405 | + # Initialise empty observation container | ||
| 406 | + obs = gammalib.GObservations() | ||
| 407 | + | ||
| 408 | + # Initialise first time and identifier | ||
| 409 | + obs_start = tstart | ||
| 410 | + obs_id = 1 | ||
| 411 | + | ||
| 412 | + # Loop over observation definition list | ||
| 413 | + for obsdef in obsdeflist: | ||
| 414 | + | ||
| 415 | + # Set pointing direction | ||
| 416 | + pntdir = gammalib.GSkyDir() | ||
| 417 | + pntdir.radec_deg(obsdef['ra'], obsdef['dec']) | ||
| 418 | + | ||
| 419 | + # Set duration (use default if not found in definition list) | ||
| 420 | + if 'duration' in obsdef: | ||
| 421 | + obs_duration = obsdef['duration'] | ||
| 422 | + else: | ||
| 423 | + obs_duration = duration | ||
| 424 | + | ||
| 425 | + # Set emin (use default if not found in definition list) | ||
| 426 | + if 'emin' in obsdef: | ||
| 427 | + obs_emin = obsdef['emin'] | ||
| 428 | + else: | ||
| 429 | + obs_emin = emin | ||
| 430 | + | ||
| 431 | + # Set emax (use default if not found in definition list) | ||
| 432 | + if 'emax' in obsdef: | ||
| 433 | + obs_emax = obsdef['emax'] | ||
| 434 | + else: | ||
| 435 | + obs_emax = emax | ||
| 436 | + | ||
| 437 | + # Set radius (use default if not found in definition list) | ||
| 438 | + if 'rad' in obsdef: | ||
| 439 | + obs_rad = obsdef['rad'] | ||
| 440 | + else: | ||
| 441 | + obs_rad = rad | ||
| 442 | + | ||
| 443 | + # Set deadc (use default if not found in definition list) | ||
| 444 | + if 'deadc' in obsdef: | ||
| 445 | + obs_deadc = obsdef['deadc'] | ||
| 446 | + else: | ||
| 447 | + obs_deadc = deadc | ||
| 448 | + | ||
| 449 | + # Set caldb (use default if not found in definition list) | ||
| 450 | + if 'caldb' in obsdef: | ||
| 451 | + obs_caldb = obsdef['caldb'] | ||
| 452 | + else: | ||
| 453 | + obs_caldb = caldb | ||
| 454 | + | ||
| 455 | + # Set irf (use default if not found in definition list) | ||
| 456 | + if 'irf' in obsdef: | ||
| 457 | + obs_irf = obsdef['irf'] | ||
| 458 | + else: | ||
| 459 | + obs_irf = irf | ||
| 460 | + | ||
| 461 | + # Generate identifier string | ||
| 462 | + id = "%6.6d" % obs_id | ||
| 463 | + | ||
| 464 | + # Set CTA observation | ||
| 465 | + obs_cta = set_obs(pntdir, tstart=obs_start, duration=obs_duration, \ | ||
| 466 | + deadc=obs_deadc, emin=obs_emin, emax=obs_emax, \ | ||
| 467 | + rad=obs_rad, irf=obs_irf, caldb=obs_caldb, id=id) | ||
| 468 | + | ||
| 469 | + # Append to container | ||
| 470 | + obs.append(obs_cta) | ||
| 471 | + | ||
| 472 | + # Update start time and identifier | ||
| 473 | + obs_start += obs_duration | ||
| 474 | + obs_id += 1 | ||
| 475 | + | ||
| 476 | + # Return observation container | ||
| 477 | + return obs | ||
| 478 | + | ||
| 479 | + | ||
| 480 | +# ============================ # | ||
| 481 | +# Set CTA observation patterns # | ||
| 482 | +# ============================ # | ||
| 483 | +def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5): | ||
| 484 | + """ | ||
| 485 | + Sets a number of standard patterns. | ||
| 486 | + | ||
| 487 | + Parameters: | ||
| 488 | + pattern - Observation pattern. Possible options are: | ||
| 489 | + - "single": single pointing | ||
| 490 | + - "four": four pointings 'offset' around pattern centre | ||
| 491 | + | ||
| 492 | + Keywords: | ||
| 493 | + ra - Right Ascension of pattern centre [deg] (default: 83.6331) | ||
| 494 | + dec - Declination of pattern centre [deg] (default: 22.0145) | ||
| 495 | + offset - Offset from pattern centre [deg] (default: 1.5) | ||
| 496 | + """ | ||
| 497 | + # Initialise observation definition list | ||
| 498 | + obsdeflist = [] | ||
| 499 | + | ||
| 500 | + # Add patterns | ||
| 501 | + if pattern == "single": | ||
| 502 | + obsdef = {'ra': ra, 'dec': dec} | ||
| 503 | + obsdeflist.append(obsdef) | ||
| 504 | + elif pattern == "four": | ||
| 505 | + # Set pattern centre | ||
| 506 | + centre = gammalib.GSkyDir() | ||
| 507 | + centre.radec_deg(ra, dec) | ||
| 508 | + | ||
| 509 | + # Append pointings | ||
| 510 | + for phi in [0.0, 90.0, 180.0, 270.0]: | ||
| 511 | + pntdir = centre.copy() | ||
| 512 | + pntdir.rotate_deg(phi, offset) | ||
| 513 | + obsdeflist.append({'ra': pntdir.ra_deg(), 'dec': pntdir.dec_deg()}) | ||
| 514 | + else: | ||
| 515 | + print("Warning: Observation pattern '"+str(pattern)+"' not recognized.") | ||
| 516 | + | ||
| 517 | + # Return observation definition list | ||
| 518 | + return obsdeflist |
pyext/setup.py.in
| @@ -121,8 +121,8 @@ setup(name='@PACKAGE_NAME@', | @@ -121,8 +121,8 @@ setup(name='@PACKAGE_NAME@', | ||
| 121 | author='Juergen Knoedlseder', | 121 | author='Juergen Knoedlseder', |
| 122 | author_email='@PACKAGE_BUGREPORT@', | 122 | author_email='@PACKAGE_BUGREPORT@', |
| 123 | url='http://gammalib.sourceforge.net/', | 123 | url='http://gammalib.sourceforge.net/', |
| 124 | - py_modules=['ctools'], | ||
| 125 | - ext_modules=[Extension('_ctools', ['ctools_wrap.cpp'], | 124 | + py_modules=['ctools.ctools', 'ctools.obsutils'], |
| 125 | + ext_modules=[Extension('ctools._ctools', ['ctools/ctools_wrap.cpp'], | ||
| 126 | include_dirs=include_dirs, | 126 | include_dirs=include_dirs, |
| 127 | library_dirs=library_dirs, | 127 | library_dirs=library_dirs, |
| 128 | libraries=libraries, | 128 | libraries=libraries, |
scripts/Makefile.am
| @@ -5,8 +5,8 @@ binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py | @@ -5,8 +5,8 @@ binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py | ||
| 5 | 5 | ||
| 6 | # Scripts to be distributed in pythondir | 6 | # Scripts to be distributed in pythondir |
| 7 | pythonscripts = cspull.py cssens.py cstsdist.py cscaldb.py \ | 7 | pythonscripts = cspull.py cssens.py cstsdist.py cscaldb.py \ |
| 8 | - obsutils.py \ | ||
| 9 | - lsmodel.py | 8 | + lsmodel.py \ |
| 9 | + obsutils.py | ||
| 10 | 10 | ||
| 11 | # Scripts for which symbolic links will be created in bindir | 11 | # Scripts for which symbolic links will be created in bindir |
| 12 | linkscripts = cspull cssens cstsdist cscaldb \ | 12 | linkscripts = cspull cssens cstsdist cscaldb \ |
scripts/cssens.py
| @@ -35,7 +35,6 @@ | @@ -35,7 +35,6 @@ | ||
| 35 | # ========================================================================== | 35 | # ========================================================================== |
| 36 | from ctools import * | 36 | from ctools import * |
| 37 | from gammalib import * | 37 | from gammalib import * |
| 38 | -import obsutils | ||
| 39 | import sys | 38 | import sys |
| 40 | import csv | 39 | import csv |
| 41 | import math | 40 | import math |
scripts/cstsdist.py
| @@ -21,7 +21,6 @@ | @@ -21,7 +21,6 @@ | ||
| 21 | # ========================================================================== | 21 | # ========================================================================== |
| 22 | from ctools import * | 22 | from ctools import * |
| 23 | from gammalib import * | 23 | from gammalib import * |
| 24 | -import obsutils | ||
| 25 | import sys | 24 | import sys |
| 26 | import csv | 25 | import csv |
| 27 | import math | 26 | import math |
scripts/obsutils.py
100755 → 100644
| 1 | -# ========================================================================== | ||
| 2 | -# This script provides a number of functions that are useful for handling | ||
| 3 | -# CTA observations. | ||
| 4 | -# | ||
| 5 | -# Copyright (C) 2011-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 | 1 | +print "The obsutils are moved under the ctools module." |
| 2 | +print "Please update your script to import from ctools, i.e:" | ||
| 3 | +print "from ctools import obsutils" | ||
| 23 | 4 | ||
| 24 | - | ||
| 25 | -# ===================== # | ||
| 26 | -# Simulate observations # | ||
| 27 | -# ===================== # | ||
| 28 | -def sim(obs, log=False, debug=False, edisp=False, seed=0, nbins=0, | ||
| 29 | - binsz=0.05, npix=200, proj="TAN", coord="GAL"): | ||
| 30 | - """ | ||
| 31 | - Simulate events for all observations in the container. | ||
| 32 | - | ||
| 33 | - Parameters: | ||
| 34 | - obs - Observation container | ||
| 35 | - Keywords: | ||
| 36 | - log - Create log file(s) | ||
| 37 | - debug - Create console dump? | ||
| 38 | - edisp - Apply energy dispersion? | ||
| 39 | - seed - Seed value for simulations (default: 0) | ||
| 40 | - nbins - Number of energy bins (default: 0=unbinned) | ||
| 41 | - binsz - Pixel size for binned simulation (deg/pixel) | ||
| 42 | - npix - Number of pixels in X and Y for binned simulation | ||
| 43 | - """ | ||
| 44 | - # Allocate ctobssim application and set parameters | ||
| 45 | - sim = ctools.ctobssim(obs) | ||
| 46 | - sim["seed"].integer(seed) | ||
| 47 | - | ||
| 48 | - # Optionally open the log file | ||
| 49 | - if log: | ||
| 50 | - sim.logFileOpen() | ||
| 51 | - | ||
| 52 | - # Optionally switch-on debugging model | ||
| 53 | - if debug: | ||
| 54 | - sim["debug"].boolean(True) | ||
| 55 | - | ||
| 56 | - # Run ctobssim application. This will loop over all observations in the | ||
| 57 | - # container and simulation the events for each observation. Note that | ||
| 58 | - # events are not added together, they still apply to each observation | ||
| 59 | - # separately. | ||
| 60 | - sim.run() | ||
| 61 | - | ||
| 62 | - # Binned option? | ||
| 63 | - if nbins > 0: | ||
| 64 | - | ||
| 65 | - # Determine common energy boundaries for observations | ||
| 66 | - emin = None | ||
| 67 | - emax = None | ||
| 68 | - for run in sim.obs(): | ||
| 69 | - run_emin = run.events().ebounds().emin().TeV() | ||
| 70 | - run_emax = run.events().ebounds().emax().TeV() | ||
| 71 | - if emin == None: | ||
| 72 | - emin = run_emin | ||
| 73 | - elif run_emin > emin: | ||
| 74 | - emin = run_emin | ||
| 75 | - if emax == None: | ||
| 76 | - emax = run_emax | ||
| 77 | - elif run_emax > emax: | ||
| 78 | - emax = run_emax | ||
| 79 | - | ||
| 80 | - # Allocate ctbin application and set parameters | ||
| 81 | - bin = ctools.ctbin(sim.obs()) | ||
| 82 | - bin["ebinalg"].string("LOG") | ||
| 83 | - bin["emin"].real(emin) | ||
| 84 | - bin["emax"].real(emax) | ||
| 85 | - bin["enumbins"].integer(nbins) | ||
| 86 | - bin["usepnt"].boolean(True) # Use pointing for map centre | ||
| 87 | - bin["nxpix"].integer(npix) | ||
| 88 | - bin["nypix"].integer(npix) | ||
| 89 | - bin["binsz"].real(binsz) | ||
| 90 | - bin["coordsys"].string(coord) | ||
| 91 | - bin["proj"].string(proj) | ||
| 92 | - | ||
| 93 | - # Optionally open the log file | ||
| 94 | - if log: | ||
| 95 | - bin.logFileOpen() | ||
| 96 | - | ||
| 97 | - # Optionally switch-on debugging model | ||
| 98 | - if debug: | ||
| 99 | - bin["debug"].boolean(True) | ||
| 100 | - | ||
| 101 | - # Run ctbin application. This will loop over all observations in | ||
| 102 | - # the container and bin the events in counts maps | ||
| 103 | - bin.run() | ||
| 104 | - | ||
| 105 | - # Make a deep copy of the observation that will be returned | ||
| 106 | - # (the ctbin object will go out of scope one the function is | ||
| 107 | - # left) | ||
| 108 | - obs = bin.obs().copy() | ||
| 109 | - | ||
| 110 | - else: | ||
| 111 | - | ||
| 112 | - # Make a deep copy of the observation that will be returned | ||
| 113 | - # (the ctobssim object will go out of scope one the function is | ||
| 114 | - # left) | ||
| 115 | - obs = sim.obs().copy() | ||
| 116 | - | ||
| 117 | - # Delete the simulation | ||
| 118 | - del sim | ||
| 119 | - | ||
| 120 | - # Return observation container | ||
| 121 | - return obs | ||
| 122 | - | ||
| 123 | - | ||
| 124 | -# ================ # | ||
| 125 | -# Fit observations # | ||
| 126 | -# ================ # | ||
| 127 | -def fit(obs, log=False, debug=False, edisp=False): | ||
| 128 | - """ | ||
| 129 | - Perform maximum likelihood fitting of observations in the container. | ||
| 130 | - | ||
| 131 | - Parameters: | ||
| 132 | - obs - Observation container | ||
| 133 | - Keywords: | ||
| 134 | - log - Create log file(s) | ||
| 135 | - debug - Create screen dump | ||
| 136 | - edisp - Apply energy dispersion? | ||
| 137 | - """ | ||
| 138 | - # Allocate ctlike application | ||
| 139 | - like = ctools.ctlike(obs) | ||
| 140 | - | ||
| 141 | - # Optionally open the log file | ||
| 142 | - if log: | ||
| 143 | - like.logFileOpen() | ||
| 144 | - | ||
| 145 | - # Optionally switch-on debugging model | ||
| 146 | - if debug: | ||
| 147 | - like["debug"].boolean(True) | ||
| 148 | - | ||
| 149 | - # Optionally apply energy dispersion | ||
| 150 | - like["edisp"].boolean(edisp) | ||
| 151 | - | ||
| 152 | - # Run ctlike application. | ||
| 153 | - like.run() | ||
| 154 | - | ||
| 155 | - # Return observations | ||
| 156 | - return like | ||
| 157 | - | ||
| 158 | - | ||
| 159 | -# ================= # | ||
| 160 | -# Create counts map # | ||
| 161 | -# ================= # | ||
| 162 | -def cntmap(obs, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ | ||
| 163 | - binsz=0.05, nxpix=200, nypix=200, \ | ||
| 164 | - outname="cntmap.fits"): | ||
| 165 | - """ | ||
| 166 | - Creates a counts map by combining the events of all observations. | ||
| 167 | - The counts map will be a summed map over all energies. | ||
| 168 | - | ||
| 169 | - Parameters: | ||
| 170 | - obs - Observation container | ||
| 171 | - Keywords: | ||
| 172 | - proj - Projection type (e.g. TAN, CAR, STG, ...) (default: TAN) | ||
| 173 | - coord - Coordinate type (GAL, CEL) (default: GAL) | ||
| 174 | - xval - Reference longitude value [deg] (default: 0.0) | ||
| 175 | - yval - Reference latitude value [deg] (default: 0.0) | ||
| 176 | - binsz - Pixel size [deg/pixel] (default: 0.05) | ||
| 177 | - nxpix - Number of pixels in X direction (default: 200) | ||
| 178 | - nypix - Number of pixels in Y direction (default: 200) | ||
| 179 | - outname - Counts map FITS filename (default: cntmap.fits) | ||
| 180 | - """ | ||
| 181 | - # Allocate counts map | ||
| 182 | - map = gammalib.GSkymap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1) | ||
| 183 | - | ||
| 184 | - # Set maximum pixel number | ||
| 185 | - maxpixel = nxpix * nypix | ||
| 186 | - | ||
| 187 | - # Fill all observations | ||
| 188 | - for run in obs: | ||
| 189 | - | ||
| 190 | - # Loop over all events | ||
| 191 | - for event in run.events(): | ||
| 192 | - | ||
| 193 | - # Determine sky pixel | ||
| 194 | - skydir = gammalib.GCTAInstDir(event.dir()).dir() | ||
| 195 | - pixel = map.dir2inx(skydir) | ||
| 196 | - | ||
| 197 | - # Set pixel | ||
| 198 | - if pixel < maxpixel: | ||
| 199 | - map[pixel] += event.counts() | ||
| 200 | - | ||
| 201 | - # Save sky map. The clobber flag is set to True, so any existing FITS | ||
| 202 | - # file will be overwritten. | ||
| 203 | - map.save(outname, True) | ||
| 204 | - | ||
| 205 | - # Return counts map | ||
| 206 | - return map | ||
| 207 | - | ||
| 208 | - | ||
| 209 | -# ================ # | ||
| 210 | -# Create model map # | ||
| 211 | -# ================ # | ||
| 212 | -def modmap(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ | ||
| 213 | - binsz=0.05, nxpix=200, nypix=200, \ | ||
| 214 | - outname="modmap.fits"): | ||
| 215 | - """ | ||
| 216 | - Make model map for a given reference energy by combining all observations. | ||
| 217 | - The model map will be evaluated for a given reference energy 'eref' and will | ||
| 218 | - be given in units of [counts/(sr MeV s)]. | ||
| 219 | - | ||
| 220 | - Parameters: | ||
| 221 | - obs - Observation container | ||
| 222 | - Keywords: | ||
| 223 | - eref - Reference energy for which model is created [TeV] (default: 0.1) | ||
| 224 | - proj - Projection type (e.g. TAN, CAR, STG, ...) (default: TAN) | ||
| 225 | - coord - Coordinate type (GAL, CEL) (default: GAL) | ||
| 226 | - xval - Reference longitude value [deg] (default: 0.0) | ||
| 227 | - yval - Reference latitude value [deg] (default: 0.0) | ||
| 228 | - binsz - Pixel size [deg/pixel] (default: 0.05) | ||
| 229 | - nxpix - Number of pixels in X direction (default: 200) | ||
| 230 | - nypix - Number of pixels in Y direction (default: 200) | ||
| 231 | - outname - Model map FITS filename (default: modmap.fits) | ||
| 232 | - """ | ||
| 233 | - # Allocate model map | ||
| 234 | - map = gammalib.GSkymap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1) | ||
| 235 | - | ||
| 236 | - # Set reference energy, time and direction. The time is not initialised and is | ||
| 237 | - # in fact not used (as the IRF is assumed to be time independent for now). | ||
| 238 | - # The sky direction is set later using the pixel values. | ||
| 239 | - energy = gammalib.GEnergy() | ||
| 240 | - time = gammalib.GTime() | ||
| 241 | - instdir = gammalib.GCTAInstDir() | ||
| 242 | - energy.TeV(eref) | ||
| 243 | - | ||
| 244 | - # Loop over all map pixels | ||
| 245 | - for pixel in range(map.npix()): | ||
| 246 | - | ||
| 247 | - # Get sky direction | ||
| 248 | - skydir = map.inx2dir(pixel) | ||
| 249 | - instdir.dir(skydir) | ||
| 250 | - | ||
| 251 | - # Create event atom for map pixel | ||
| 252 | - atom = gammalib.GCTAEventAtom() | ||
| 253 | - atom.dir(instdir) | ||
| 254 | - atom.energy(energy) | ||
| 255 | - atom.time(time) | ||
| 256 | - | ||
| 257 | - # Initialise model value | ||
| 258 | - value = 0.0 | ||
| 259 | - | ||
| 260 | - # Loop over all observations | ||
| 261 | - for run in obs: | ||
| 262 | - value += obs.models().eval(atom, run) | ||
| 263 | - | ||
| 264 | - # Set map value | ||
| 265 | - map[pixel] = value | ||
| 266 | - | ||
| 267 | - # Save sky map | ||
| 268 | - map.save(outname, True) | ||
| 269 | - | ||
| 270 | - # Return model map | ||
| 271 | - return map | ||
| 272 | - | ||
| 273 | - | ||
| 274 | -# ===================================== # | ||
| 275 | -# Set one CTA observation (old version) # | ||
| 276 | -# ===================================== # | ||
| 277 | -def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ | ||
| 278 | - emin=0.1, emax=100.0, rad=5.0, \ | ||
| 279 | - irf="cta_dummy_irf", caldb="dummy"): | ||
| 280 | - """ | ||
| 281 | - Obsolete function, use set_obs instead. | ||
| 282 | - """ | ||
| 283 | - # Print warning | ||
| 284 | - print("Warning: obsutils.set is obsolete, use obsutils.set_obs instead.") | ||
| 285 | - | ||
| 286 | - # Call new function | ||
| 287 | - obs_cta = set_obs(pntdir, tstart=tstart, duration=duration, deadc=deadc, \ | ||
| 288 | - emin=emin, emax=emax, rad=rad, \ | ||
| 289 | - irf=irf, caldb=caldb) | ||
| 290 | - | ||
| 291 | - # Return CTA observation | ||
| 292 | - return obs_cta | ||
| 293 | - | ||
| 294 | - | ||
| 295 | -# ======================= # | ||
| 296 | -# Set one CTA observation # | ||
| 297 | -# ======================= # | ||
| 298 | -def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ | ||
| 299 | - emin=0.1, emax=100.0, rad=5.0, \ | ||
| 300 | - irf="cta_dummy_irf", caldb="dummy", id="000000"): | ||
| 301 | - """ | ||
| 302 | - Returns a single CTA observation containing an empty CTA event list. | ||
| 303 | - By looping over this function you can add CTA observations to the | ||
| 304 | - observation container. | ||
| 305 | - | ||
| 306 | - Parameters: | ||
| 307 | - pntdir - Pointing direction [GSkyDir] | ||
| 308 | - Keywords: | ||
| 309 | - tstart - Start time [seconds] (default: 0.0) | ||
| 310 | - duration - Duration of observation [seconds] (default: 1800.0) | ||
| 311 | - deadc - Deadtime correction factor (default: 0.95) | ||
| 312 | - emin - Minimum event energy [TeV] (default: 0.1) | ||
| 313 | - emax - Maximum event energy [TeV] (default: 100.0) | ||
| 314 | - rad - ROI radius used for analysis [deg] (default: 5.0) | ||
| 315 | - irf - Instrument response function (default: cta_dummy_irf) | ||
| 316 | - caldb - Calibration database path (default: "dummy") | ||
| 317 | - """ | ||
| 318 | - # Allocate CTA observation | ||
| 319 | - obs_cta = gammalib.GCTAObservation() | ||
| 320 | - | ||
| 321 | - # Set calibration database | ||
| 322 | - db = gammalib.GCaldb() | ||
| 323 | - if (gammalib.dir_exists(caldb)): | ||
| 324 | - db.rootdir(caldb) | ||
| 325 | - else: | ||
| 326 | - db.open("cta", caldb) | ||
| 327 | - | ||
| 328 | - # Set pointing direction | ||
| 329 | - pnt = gammalib.GCTAPointing() | ||
| 330 | - pnt.dir(pntdir) | ||
| 331 | - obs_cta.pointing(pnt) | ||
| 332 | - | ||
| 333 | - # Set ROI | ||
| 334 | - roi = gammalib.GCTARoi() | ||
| 335 | - instdir = gammalib.GCTAInstDir() | ||
| 336 | - instdir.dir(pntdir) | ||
| 337 | - roi.centre(instdir) | ||
| 338 | - roi.radius(rad) | ||
| 339 | - | ||
| 340 | - # Set GTI | ||
| 341 | - gti = gammalib.GGti() | ||
| 342 | - gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration)) | ||
| 343 | - | ||
| 344 | - # Set energy boundaries | ||
| 345 | - ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \ | ||
| 346 | - gammalib.GEnergy(emax, "TeV")) | ||
| 347 | - | ||
| 348 | - # Allocate event list | ||
| 349 | - events = gammalib.GCTAEventList() | ||
| 350 | - events.roi(roi) | ||
| 351 | - events.gti(gti) | ||
| 352 | - events.ebounds(ebounds) | ||
| 353 | - obs_cta.events(events) | ||
| 354 | - | ||
| 355 | - # Set instrument response | ||
| 356 | - obs_cta.response(irf, db) | ||
| 357 | - | ||
| 358 | - # Set ontime, livetime, and deadtime correction factor | ||
| 359 | - obs_cta.ontime(duration) | ||
| 360 | - obs_cta.livetime(duration*deadc) | ||
| 361 | - obs_cta.deadc(deadc) | ||
| 362 | - obs_cta.id(id) | ||
| 363 | - | ||
| 364 | - # Return CTA observation | ||
| 365 | - return obs_cta | ||
| 366 | - | ||
| 367 | - | ||
| 368 | -# ============================ # | ||
| 369 | -# Set list of CTA observations # | ||
| 370 | -# ============================ # | ||
| 371 | -def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \ | ||
| 372 | - emin=0.1, emax=100.0, rad=5.0, \ | ||
| 373 | - irf="cta_dummy_irf", caldb="dummy"): | ||
| 374 | - """ | ||
| 375 | - Returns an observation container filled with a list of CTA observations. | ||
| 376 | - The list is defined by the obsdeflist parameter which is a dictionnary | ||
| 377 | - containing the mandatory keywords 'ra' and 'dec' that specify the | ||
| 378 | - pointing direction for a given observation. Optional keyword give control | ||
| 379 | - over other observation proposerties, such as duration, deadtime correction, | ||
| 380 | - energy range, etc. If an optional keyword is not specified, the function | ||
| 381 | - keyword is used instead. | ||
| 382 | - | ||
| 383 | - Parameters: | ||
| 384 | - obsdeflist - Observation definition list [{'ra': x.xx, 'dec': x.xx}] | ||
| 385 | - The directory can take the following optional keywords: | ||
| 386 | - - duration: Duration of CTA observation [seconds] | ||
| 387 | - - deadc: Deadtime correction factor | ||
| 388 | - - emin: Minimum event energy [TeV] | ||
| 389 | - - emax: Maximum event energy [TeV] | ||
| 390 | - - rad: ROI radius used for analysis [deg] | ||
| 391 | - - irf: Instrument response function | ||
| 392 | - - caldb: Calibration database path | ||
| 393 | - Optional keywords overwrite keywords specified in the | ||
| 394 | - function call. | ||
| 395 | - Keywords: | ||
| 396 | - tstart - Start time [seconds] (default: 0.0) | ||
| 397 | - duration - Duration of one CTA observation [seconds] (default: 1800.0) | ||
| 398 | - deadc - Deadtime correction factor (default: 0.95) | ||
| 399 | - emin - Minimum event energy [TeV] (default: 0.1) | ||
| 400 | - emax - Maximum event energy [TeV] (default: 100.0) | ||
| 401 | - rad - ROI radius used for analysis [deg] (default: 5.0) | ||
| 402 | - irf - Instrument response function (default: cta_dummy_irf) | ||
| 403 | - caldb - Calibration database path (default: "dummy") | ||
| 404 | - """ | ||
| 405 | - # Initialise empty observation container | ||
| 406 | - obs = gammalib.GObservations() | ||
| 407 | - | ||
| 408 | - # Initialise first time and identifier | ||
| 409 | - obs_start = tstart | ||
| 410 | - obs_id = 1 | ||
| 411 | - | ||
| 412 | - # Loop over observation definition list | ||
| 413 | - for obsdef in obsdeflist: | ||
| 414 | - | ||
| 415 | - # Set pointing direction | ||
| 416 | - pntdir = gammalib.GSkyDir() | ||
| 417 | - pntdir.radec_deg(obsdef['ra'], obsdef['dec']) | ||
| 418 | - | ||
| 419 | - # Set duration (use default if not found in definition list) | ||
| 420 | - if 'duration' in obsdef: | ||
| 421 | - obs_duration = obsdef['duration'] | ||
| 422 | - else: | ||
| 423 | - obs_duration = duration | ||
| 424 | - | ||
| 425 | - # Set emin (use default if not found in definition list) | ||
| 426 | - if 'emin' in obsdef: | ||
| 427 | - obs_emin = obsdef['emin'] | ||
| 428 | - else: | ||
| 429 | - obs_emin = emin | ||
| 430 | - | ||
| 431 | - # Set emax (use default if not found in definition list) | ||
| 432 | - if 'emax' in obsdef: | ||
| 433 | - obs_emax = obsdef['emax'] | ||
| 434 | - else: | ||
| 435 | - obs_emax = emax | ||
| 436 | - | ||
| 437 | - # Set radius (use default if not found in definition list) | ||
| 438 | - if 'rad' in obsdef: | ||
| 439 | - obs_rad = obsdef['rad'] | ||
| 440 | - else: | ||
| 441 | - obs_rad = rad | ||
| 442 | - | ||
| 443 | - # Set deadc (use default if not found in definition list) | ||
| 444 | - if 'deadc' in obsdef: | ||
| 445 | - obs_deadc = obsdef['deadc'] | ||
| 446 | - else: | ||
| 447 | - obs_deadc = deadc | ||
| 448 | - | ||
| 449 | - # Set caldb (use default if not found in definition list) | ||
| 450 | - if 'caldb' in obsdef: | ||
| 451 | - obs_caldb = obsdef['caldb'] | ||
| 452 | - else: | ||
| 453 | - obs_caldb = caldb | ||
| 454 | - | ||
| 455 | - # Set irf (use default if not found in definition list) | ||
| 456 | - if 'irf' in obsdef: | ||
| 457 | - obs_irf = obsdef['irf'] | ||
| 458 | - else: | ||
| 459 | - obs_irf = irf | ||
| 460 | - | ||
| 461 | - # Generate identifier string | ||
| 462 | - id = "%6.6d" % obs_id | ||
| 463 | - | ||
| 464 | - # Set CTA observation | ||
| 465 | - obs_cta = set_obs(pntdir, tstart=obs_start, duration=obs_duration, \ | ||
| 466 | - deadc=obs_deadc, emin=obs_emin, emax=obs_emax, \ | ||
| 467 | - rad=obs_rad, irf=obs_irf, caldb=obs_caldb, id=id) | ||
| 468 | - | ||
| 469 | - # Append to container | ||
| 470 | - obs.append(obs_cta) | ||
| 471 | - | ||
| 472 | - # Update start time and identifier | ||
| 473 | - obs_start += obs_duration | ||
| 474 | - obs_id += 1 | ||
| 475 | - | ||
| 476 | - # Return observation container | ||
| 477 | - return obs | ||
| 478 | - | ||
| 479 | - | ||
| 480 | -# ============================ # | ||
| 481 | -# Set CTA observation patterns # | ||
| 482 | -# ============================ # | ||
| 483 | -def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5): | ||
| 484 | - """ | ||
| 485 | - Sets a number of standard patterns. | ||
| 486 | - | ||
| 487 | - Parameters: | ||
| 488 | - pattern - Observation pattern. Possible options are: | ||
| 489 | - - "single": single pointing | ||
| 490 | - - "four": four pointings 'offset' around pattern centre | ||
| 491 | - | ||
| 492 | - Keywords: | ||
| 493 | - ra - Right Ascension of pattern centre [deg] (default: 83.6331) | ||
| 494 | - dec - Declination of pattern centre [deg] (default: 22.0145) | ||
| 495 | - offset - Offset from pattern centre [deg] (default: 1.5) | ||
| 496 | - """ | ||
| 497 | - # Initialise observation definition list | ||
| 498 | - obsdeflist = [] | ||
| 499 | - | ||
| 500 | - # Add patterns | ||
| 501 | - if pattern == "single": | ||
| 502 | - obsdef = {'ra': ra, 'dec': dec} | ||
| 503 | - obsdeflist.append(obsdef) | ||
| 504 | - elif pattern == "four": | ||
| 505 | - # Set pattern centre | ||
| 506 | - centre = gammalib.GSkyDir() | ||
| 507 | - centre.radec_deg(ra, dec) | ||
| 508 | - | ||
| 509 | - # Append pointings | ||
| 510 | - for phi in [0.0, 90.0, 180.0, 270.0]: | ||
| 511 | - pntdir = centre.copy() | ||
| 512 | - pntdir.rotate_deg(phi, offset) | ||
| 513 | - obsdeflist.append({'ra': pntdir.ra_deg(), 'dec': pntdir.dec_deg()}) | ||
| 514 | - else: | ||
| 515 | - print("Warning: Observation pattern '"+str(pattern)+"' not recognized.") | ||
| 516 | - | ||
| 517 | - # Return observation definition list | ||
| 518 | - return obsdeflist | 5 | +from ctools.obsutils import * |
test/Makefile.am
| @@ -8,7 +8,7 @@ MACOSX_DEPLOYMENT_TARGET=@CONFIGURE_MACOSX_DEPLOYMENT_TARGET@ | @@ -8,7 +8,7 @@ MACOSX_DEPLOYMENT_TARGET=@CONFIGURE_MACOSX_DEPLOYMENT_TARGET@ | ||
| 8 | if WITH_PYTHON | 8 | if WITH_PYTHON |
| 9 | TEST_CSCRIPTS = test_cscripts.sh | 9 | TEST_CSCRIPTS = test_cscripts.sh |
| 10 | TEST_PYTHON_SCRIPT = test_python.py | 10 | TEST_PYTHON_SCRIPT = test_python.py |
| 11 | - TEST_PYTHON_ENV = PYTHONPATH=$(top_builddir)/pyext:$(top_builddir)/pyext/build/ctools:$(top_builddir)/scripts:$(pythondir):$(pyexecdir):$(PYTHONPATH) | 11 | + TEST_PYTHON_ENV = PYTHONPATH=$(top_builddir)/pyext:$(top_builddir)/pyext/build/ctools/ctools:$(top_builddir)/scripts:$(pythondir):$(pyexecdir):$(PYTHONPATH) |
| 12 | endif | 12 | endif |
| 13 | 13 | ||
| 14 | # Set test environment | 14 | # Set test environment |