Commit ad7ffe86388e9601625556421f1cda53a6ce813d
Merge branch 'integration' into 1244-extendgtlike
Conflicts: scripts/obsutils.py
Showing
67 changed files
with
1858 additions
and
1947 deletions
.gitignore
| ... | ... | @@ -39,19 +39,20 @@ depcomp |
| 39 | 39 | test-driver |
| 40 | 40 | testconf.sh |
| 41 | 41 | |
| 42 | -doc/html/doxygen | |
| 42 | +doc/html | |
| 43 | 43 | doc/man |
| 44 | 44 | src/ctools-setup |
| 45 | 45 | pyext/build/ |
| 46 | 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 | 49 | src/ctbin/ctbin |
| 50 | 50 | src/ctlike/ctlike |
| 51 | 51 | src/ctmodel/ctmodel |
| 52 | 52 | src/ctobssim/ctobssim |
| 53 | 53 | src/ctselect/ctselect |
| 54 | 54 | src/ctskymap/ctskymap |
| 55 | +src/ctexpcube/ctexpcube | |
| 55 | 56 | test/pfiles/ |
| 56 | 57 | test/*.fits |
| 57 | 58 | test/*.dat |
| ... | ... | @@ -61,20 +62,3 @@ test/test_cscripts.sh |
| 61 | 62 | examples/*.fits |
| 62 | 63 | examples/*.xml |
| 63 | 64 | examples/*.dat |
| 64 | - | |
| 65 | -# Sphinx stuff | |
| 66 | -doc/html/.buildinfo | |
| 67 | -doc/html/.doctrees/ | |
| 68 | -doc/html/_sources/ | |
| 69 | -doc/html/_static/ | |
| 70 | -doc/html/_images/ | |
| 71 | -doc/html/contents.html | |
| 72 | -doc/html/develop.html | |
| 73 | -doc/html/download.html | |
| 74 | -doc/html/genindex.html | |
| 75 | -doc/html/objects.inv | |
| 76 | -doc/html/search.html | |
| 77 | -doc/html/searchindex.js | |
| 78 | -doc/html/test.html | |
| 79 | -doc/html/index.html | |
| 80 | -doc/html/user_manual/ | ... | ... |
ChangeLog
| 1 | -2014-07-09 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu> | |
| 1 | +2014-07-10 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu> | |
| 2 | 2 | |
| 3 | 3 | * Version 0.8.0 released |
| 4 | 4 | ======================== |
| 5 | 5 | |
| 6 | + Add ctexpcube tool | |
| 6 | 7 | Support energy dispersion in ctobssim, ctlike and ctmodel |
| 7 | 8 | Correct ctobssim for thread safety |
| 8 | 9 | Enhance documentation | ... | ... |
Makefile.am
| ... | ... | @@ -25,14 +25,20 @@ doc_DATA = \ |
| 25 | 25 | README |
| 26 | 26 | docdir = ${datadir}/doc/${PACKAGE} |
| 27 | 27 | |
| 28 | -install-data-local: $(HTML_INSTALL_DATA_TARGETS) | |
| 28 | +#install-data-local: $(HTML_INSTALL_DATA_TARGETS) | |
| 29 | +#uninstall-local: $(HTML_UNINSTALL_DATA_TARGETS) | |
| 29 | 30 | |
| 30 | -uninstall-local: $(HTML_UNINSTALL_DATA_TARGETS) | |
| 31 | +install-data-local: $(DOXYGEN_INSTALL_DATA_TARGETS) sphinx-install | |
| 32 | + | |
| 33 | +uninstall-local: $(DOXYGEN_UNINSTALL_DATA_TARGETS) | |
| 31 | 34 | |
| 32 | 35 | # |
| 33 | 36 | # Additional rules |
| 34 | 37 | # |
| 35 | 38 | |
| 39 | +# Rules | |
| 40 | +.PHONY: doxygen-run doxygen-doc $(DX_PS_GOAL) $(DX_PDF_GOAL) | |
| 41 | + | |
| 36 | 42 | # Print out all rules |
| 37 | 43 | help: |
| 38 | 44 | @echo "make targets:"; \ |
| ... | ... | @@ -227,25 +233,30 @@ doxygen: doxygen-run $(DX_PS_GOAL) $(DX_PDF_GOAL) |
| 227 | 233 | |
| 228 | 234 | # Install Doxygen files |
| 229 | 235 | DOCDIRS = html/doxygen man/man3 |
| 236 | +DOXYGEN_INSTALL_DATA_TARGETS = doxygen-install | |
| 230 | 237 | doxygen-install: |
| 231 | 238 | $(mkinstalldirs) $(DESTDIR)$(docdir) |
| 232 | 239 | @for dir in $(DOCDIRS) ; do \ |
| 233 | - $(mkinstalldirs) $(DESTDIR)$(docdir)/$$dir && \ | |
| 234 | - for file in $(srcdir)/doc/$$dir/*.* ; do \ | |
| 235 | - echo "$(INSTALL_DATA) $$file $(DESTDIR)$(docdir)/$$dir" ; \ | |
| 236 | - $(INSTALL_DATA) "$$file" $(DESTDIR)$(docdir)/$$dir ; \ | |
| 237 | - done ; \ | |
| 240 | + if test -e $(srcdir)/doc/$$dir ; then \ | |
| 241 | + $(mkinstalldirs) $(DESTDIR)$(docdir)/$$dir && \ | |
| 242 | + for file in $(srcdir)/doc/$$dir/*.* ; do \ | |
| 243 | + echo "$(INSTALL_DATA) $$file $(DESTDIR)$(docdir)/$$dir" ; \ | |
| 244 | + $(INSTALL_DATA) "$$file" $(DESTDIR)$(docdir)/$$dir ; \ | |
| 245 | + done ; \ | |
| 246 | + fi ; \ | |
| 238 | 247 | done |
| 239 | 248 | |
| 240 | 249 | # Uninstall Doxygen files |
| 250 | +DOXYGEN_UNINSTALL_DATA_TARGETS = doxygen-uninstall | |
| 241 | 251 | doxygen-uninstall: |
| 242 | 252 | for dir in $(DOCDIRS) ; do \ |
| 243 | 253 | rm -f -r $(DESTDIR)$(docdir)/$$dir ; \ |
| 244 | 254 | done |
| 245 | 255 | |
| 246 | 256 | @DX_DOCDIR@/@PACKAGE@.tag: $(DX_CONFIG) $(pkginclude_HEADERS) |
| 247 | - rm -rf @DX_DOCDIR@ | |
| 248 | - $(DX_ENV) $(DX_DOXYGEN) $(srcdir)/$(DX_CONFIG) | |
| 257 | + rm -rf @DX_DOCDIR@ | |
| 258 | + $(mkinstalldirs) $(srcdir)/doc/html/doxygen | |
| 259 | + $(DX_ENV) $(DX_DOXYGEN) $(srcdir)/$(DX_CONFIG) | |
| 249 | 260 | |
| 250 | 261 | DX_CLEANFILES = \ |
| 251 | 262 | @DX_DOCDIR@/@PACKAGE@.tag \ |
| ... | ... | @@ -267,7 +278,8 @@ endif #DX_COND_doc |
| 267 | 278 | # |
| 268 | 279 | |
| 269 | 280 | # Install HTML files |
| 270 | -HTMLDIRS = html html/images | |
| 281 | +#HTMLDIRS = html | |
| 282 | +HTMLDIRS = | |
| 271 | 283 | HTML_INSTALL_DATA_TARGETS = install-data-html |
| 272 | 284 | install-data-html: |
| 273 | 285 | $(mkinstalldirs) $(DESTDIR)$(docdir) | ... | ... |
NEWS
| ... | ... | @@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-08-00 |
| 2 | 2 | |
| 3 | 3 | Author(s) : Juergen Knoedlseder |
| 4 | 4 | |
| 5 | -9 July 2014 | |
| 5 | +10 July 2014 | |
| 6 | 6 | |
| 7 | 7 | |
| 8 | 8 | 1. Introduction |
| ... | ... | @@ -65,24 +65,33 @@ the use for 2D response functions in FITS files. |
| 65 | 65 | ctmodel now also supports the handling of energy dispersion. Energy |
| 66 | 66 | dispersion can be switched on or off using the hidden edisp parameter. |
| 67 | 67 | |
| 68 | +10. ctexpcube - CTA exposure cube | |
| 69 | +--------------------------------- | |
| 70 | +The first version of the ctexpcube tool is included in the release. | |
| 71 | +ctexpcube computes the exposure cube for an observation definition | |
| 72 | +file. The cube definition may be taken from a counts map or may be | |
| 73 | +specified by user parameters. | |
| 74 | + | |
| 68 | 75 | |
| 69 | -10. Examples | |
| 76 | +11. Examples | |
| 70 | 77 | ------------ |
| 71 | 78 | None |
| 72 | 79 | |
| 73 | 80 | |
| 74 | -11. Scripts | |
| 81 | +12. Scripts | |
| 75 | 82 | ----------- |
| 76 | 83 | Add logL and event statistics to cspull output file. |
| 77 | 84 | |
| 78 | 85 | Convert all tabs into spaces. |
| 79 | 86 | |
| 87 | +The obsutils script has now been moved into the ctools module namespace. | |
| 88 | + | |
| 80 | 89 | |
| 81 | -12. caldb | |
| 90 | +13. caldb | |
| 82 | 91 | --------- |
| 83 | 92 | None |
| 84 | 93 | |
| 85 | 94 | |
| 86 | -13. models | |
| 95 | +14. models | |
| 87 | 96 | ---------- |
| 88 | 97 | None | ... | ... |
configure.ac
| ... | ... | @@ -630,7 +630,7 @@ else |
| 630 | 630 | |
| 631 | 631 | # Do we have wrappers (ctools_wrap.cpp and ctools.py)? |
| 632 | 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 | 634 | [has_wrappers="no"]) |
| 635 | 635 | if test "x$has_wrappers" = "xno"; then |
| 636 | 636 | AC_MSG_WARN([Python wrapper missing. Requires swig for wrapper generation.]) |
| ... | ... | @@ -669,7 +669,7 @@ if test "x$enable_python" = "xyes"; then |
| 669 | 669 | |
| 670 | 670 | else |
| 671 | 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 | 673 | AC_MSG_WARN([Python binding for ctools cannot be built]) |
| 674 | 674 | AC_MSG_WARN([Either use a ctools distribution with ctools_wrap.cpp and ctools.py shipped or install swig on your system]) |
| 675 | 675 | AC_MSG_WARN([swig can be obtained from http://www.swig.org/]) |
| ... | ... | @@ -882,6 +882,7 @@ AC_CONFIG_FILES([Makefile |
| 882 | 882 | src/ctobssim/Makefile |
| 883 | 883 | src/ctskymap/Makefile |
| 884 | 884 | src/ctselect/Makefile |
| 885 | + src/ctexpcube/Makefile | |
| 885 | 886 | src/ctools-setup |
| 886 | 887 | test/Makefile |
| 887 | 888 | test/test_cscripts.sh | ... | ... |
doc/Doxyfile
| ... | ... | @@ -1553,7 +1553,7 @@ TAGFILES = |
| 1553 | 1553 | # When a file name is specified after GENERATE_TAGFILE, doxygen will create |
| 1554 | 1554 | # a tag file that is based on the input files it reads. |
| 1555 | 1555 | |
| 1556 | -GENERATE_TAGFILE = | |
| 1556 | +GENERATE_TAGFILE = doc/html/doxygen/ctools.tag | |
| 1557 | 1557 | |
| 1558 | 1558 | # If the ALLEXTERNALS tag is set to YES all external classes will be listed |
| 1559 | 1559 | # in the class index. If set to NO only the inherited external classes | ... | ... |
doc/html/documentation.html deleted
| 1 | -<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"> | |
| 2 | -<html><head> | |
| 3 | -<meta content="text/html; charset=ISO-8859-1" http-equiv="content-type"><title>Documentation</title> | |
| 4 | - | |
| 5 | -<link rel="stylesheet" type="text/css" media="screen,print" href="ctools.css"> | |
| 6 | -</head> | |
| 7 | -<body> | |
| 8 | -<p class="Header">Documentation</p> | |
| 9 | -<br> | |
| 10 | -<p style="padding-top: 0pt; text-align: left;" class="Title">ctools | |
| 11 | -references | |
| 12 | -</p> | |
| 13 | -<ul class="Body"> | |
| 14 | -<li><a href="doc/ctbin.html">ctbin</a></li> | |
| 15 | -<li><a href="doc/ctlike.html">ctlike</a></li> | |
| 16 | -<li><a href="doc/ctmodel.html">ctmodel</a></li> | |
| 17 | -<li><a href="doc/ctobssim.html">ctobssim</a></li> | |
| 18 | -<li><a href="doc/ctselect.html">ctselect</a></li> | |
| 19 | -<li><a href="doc/ctskymap.html">ctskymap</a></li> | |
| 20 | -</ul> | |
| 21 | -<br> | |
| 22 | -<p style="padding-top: 0pt; text-align: left;" class="Title">cscripts | |
| 23 | -references | |
| 24 | -</p> | |
| 25 | -<ul class="Body"> | |
| 26 | -<li><a href="doc/cspull.html">cspull</a></li> | |
| 27 | -<li><a href="doc/cssens.html">cssens</a></li> | |
| 28 | -<li><a href="doc/cstsdist.html">cstsdist</a></li> | |
| 29 | -</ul> | |
| 30 | -<br> | |
| 31 | -<hr size="1"> | |
| 32 | -<p class="Body">Last updated by: Jürgen | |
| 33 | -Knödlseder (21 December 2012) | |
| 34 | -</p> | |
| 35 | -</body></html> | |
| 36 | 0 | \ No newline at end of file |
doc/html/home.html deleted
| 1 | -<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"> | |
| 2 | -<html><head> | |
| 3 | -<meta content="text/html; charset=ISO-8859-1" http-equiv="content-type"><title>Home</title> | |
| 4 | - | |
| 5 | -<link rel="stylesheet" type="text/css" media="screen,print" href="ctools.css"> | |
| 6 | -</head> | |
| 7 | -<body> | |
| 8 | -<p class="Header">ctools</p> | |
| 9 | -<p style="padding-top: 0pt; text-align: left;" class="Body">ctools | |
| 10 | -are a set of ftools-like executables needed for the scientific analysis | |
| 11 | -of Cherenkov Telescope Array data. ctools is also a Python module | |
| 12 | -that allows for interactive data analysis and building of analysis | |
| 13 | -scripts and pipelines. ctools includes also an observation simulator | |
| 14 | -to enable the scientific simulation of future CTA observations.</p> | |
| 15 | -<br> | |
| 16 | -<p style="padding-top: 0pt; text-align: left;" class="Body">ctools | |
| 17 | -are based on <a href="http://gammalib.sourceforge.net/" target="_blank">GammaLib</a>, a versatile toolbox for | |
| 18 | -the high-level analysis of astronomical gamma-ray data. Besides CTA, | |
| 19 | -GammaLib supports also the analysis of Fermi-LAT data, and extensions | |
| 20 | -to support further gamma-ray instruments are planned. This enables a | |
| 21 | -simultaneous and coherent multi-instrument analysis of high-energy | |
| 22 | -sources in the Universe.</p> | |
| 23 | -<br> | |
| 24 | -<p style="padding-top: 0pt; text-align: left;" class="Body">ctools | |
| 25 | -and GammaLib is free software distributed under the <a href="http://www.gnu.org/licenses/gpl.html" target="_blank">GPL | |
| 26 | -license version 3</a>. </p> | |
| 27 | -<br> | |
| 28 | -<table style="width: 100%; text-align: left; margin-left: auto; margin-right: auto;" border="0" cellpadding="2" cellspacing="2"> | |
| 29 | -<tbody> | |
| 30 | -<tr> | |
| 31 | -<td style="text-align: center;"></td> | |
| 32 | -<td style="text-align: center; width: 600px;" class="Body"><img style="width: 600px; height: 498px;" alt="Example simulation" src="images/survey-pop2.jpg"><br> | |
| 33 | -Simulation | |
| 34 | -of galactic plane scan (Dubus et al. 2012).</td> | |
| 35 | -<td style="text-align: center;"></td> | |
| 36 | -</tr> | |
| 37 | -</tbody> | |
| 38 | -</table> | |
| 39 | -<br> | |
| 40 | -<hr size="1"> | |
| 41 | -<p class="Body">Last updated by: Jürgen | |
| 42 | -Knödlseder (4 April 2012)<br> | |
| 43 | -</p> | |
| 44 | -</body></html> | |
| 45 | 0 | \ No newline at end of file |
doc/html/images/cntmap-fits.jpg deleted
47 KB
doc/html/images/cntmap-map.jpg deleted
22.7 KB
doc/html/images/dskeys.jpg deleted
99.6 KB
doc/html/installation.html deleted
| 1 | -<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"> | |
| 2 | -<html><head> | |
| 3 | -<meta http-equiv="content-type" content="text/html; charset=ISO-8859-1"><title>Installation</title> | |
| 4 | - | |
| 5 | -<link href="ctools.css" media="screen,print" type="text/css" rel="stylesheet"></head> | |
| 6 | -<body> | |
| 7 | -<p class="Header">Installation</p> | |
| 8 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 9 | -<a href="#Step_1:_Install_GammaLib">Step 1: Install | |
| 10 | -GammaLib</a> | |
| 11 | -</p> | |
| 12 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 13 | -<a href="#Step_2:_Setup_GammaLib_environment">Step 2: | |
| 14 | -Setup | |
| 15 | -GammaLib environment</a> | |
| 16 | -</p> | |
| 17 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 18 | -<a href="#Step_3:_Install_ctatools">Step 3: Install | |
| 19 | -ctools</a> | |
| 20 | -</p> | |
| 21 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 22 | -<a href="#Step_4:_Setup_ctatools_environment">Step 4: | |
| 23 | -Setup | |
| 24 | -ctools environment</a> | |
| 25 | -</p> | |
| 26 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><a href="#Known_Problems">Known Problems</a></p> | |
| 27 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><a href="#Getting_Help">Getting Help</a></p> | |
| 28 | -<br> | |
| 29 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><span style="font-weight: bold;"> | |
| 30 | -Note:</span> These installation instructions apply to | |
| 31 | -<span style="font-family: Courier New,Courier,monospace;">ctools-00-04-00</span> | |
| 32 | -or later. | |
| 33 | -</p> | |
| 34 | -<br> | |
| 35 | -<p class="Title"><a name="Step_1:_Install_GammaLib"></a>Step | |
| 36 | -1: Install GammaLib</p> | |
| 37 | -<p style="padding-top: 0pt; text-align: left;" class="Body">ctools | |
| 38 | -depends on GammaLib and thus GammaLib needs to be installed on your | |
| 39 | -system prior to installing ctools. If you do not have GammaLib | |
| 40 | -installed, or if you need to upgrade GammaLib, fetch the latest copy | |
| 41 | -from <a href="https://sourceforge.net/projects/gammalib/" target="_blank">https://sourceforge.net/projects/gammalib/</a> | |
| 42 | -and | |
| 43 | -follow the GammaLib installation instructions on | |
| 44 | -<a href="http://gammalib.sourceforge.net/" target="_blank">http://gammalib.sourceforge.net/</a>. | |
| 45 | -<br> | |
| 46 | -<br> | |
| 47 | -</p> | |
| 48 | -<p class="Title"><a name="Step_2:_Setup_GammaLib_environment"></a>Step | |
| 49 | -2: Setup GammaLib environment</p> | |
| 50 | -<p style="padding-top: 0pt; text-align: left;" class="Body">Before | |
| 51 | -building and installing ctools, you have to configure GammaLib by | |
| 52 | -setting up some | |
| 53 | -environment variables. | |
| 54 | -This will be done automatically by an initialisation script that is | |
| 55 | -found in the <span style="font-family: Courier New,Courier,monospace;">bin</span> | |
| 56 | -directory of the GammaLib directory (we call here "GammaLib | |
| 57 | -directory" the directory into which GammaLib has been | |
| 58 | -installed). Assuming that you have installed GammaLib in the default | |
| 59 | -directory | |
| 60 | -<span style="font-family: Courier New,Courier,monospace;">/usr/local/gamma</span> | |
| 61 | -you need to add the following to your <span style="font-family: Courier New,Courier,monospace;">$HOME/.bashrc</span> | |
| 62 | -or <span style="font-family: Courier New,Courier,monospace;">$HOME/.profile</span> | |
| 63 | -script on a Linux | |
| 64 | -machine: | |
| 65 | -</p> | |
| 66 | -<p class="Code">export GAMMALIB=/usr/local/gamma<br> | |
| 67 | -source $GAMMALIB/bin/gammalib-init.sh | |
| 68 | -</p> | |
| 69 | -<p style="padding-top: 0pt; text-align: left;" class="Body">If | |
| 70 | -you use C shell or a variant then add the following to your | |
| 71 | -<span style="font-family: Courier New,Courier,monospace;">$HOME/.cshrc</span> | |
| 72 | -or <span style="font-family: Courier New,Courier,monospace;">$HOME/.tcshrc</span> | |
| 73 | -script: | |
| 74 | -</p> | |
| 75 | -<p class="Code">setenv GAMMALIB /usr/local/gamma<br> | |
| 76 | -source $GAMMALIB/bin/gammalib-init.csh </p> | |
| 77 | -<br> | |
| 78 | -<br> | |
| 79 | -<p class="Title"><a name="Step_3:_Install_ctatools"></a>Step | |
| 80 | -3: Install ctools</p> | |
| 81 | -<p style="padding-top: 0pt; text-align: left;" class="Body">After | |
| 82 | -downloading the ctools tarball, save it in an appropriate location | |
| 83 | -(for example <span style="font-family: Courier New,Courier,monospace;">/usr/local/gamma/src</span>), | |
| 84 | -and type</p> | |
| 85 | -<p class="Code"> | |
| 86 | -> tar xvfz ctools-00-04-00.tar.gz | |
| 87 | -</p> | |
| 88 | -<p style="padding-top: 0pt; text-align: left;" class="Body">Step | |
| 89 | -in the created directory using | |
| 90 | -</p> | |
| 91 | -<p class="Code"> | |
| 92 | -> cd ctools-00-04-00 | |
| 93 | -</p> | |
| 94 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 95 | -and build and install ctools by typing</p> | |
| 96 | -<p class="Code"> | |
| 97 | -> ./configure<br> | |
| 98 | -> make<br> | |
| 99 | -> make check<br> | |
| 100 | -> make install | |
| 101 | -</p> | |
| 102 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 103 | -at the operating system prompt. | |
| 104 | -The <span style="font-family: Courier New,Courier,monospace;">./configure</span> | |
| 105 | -command customizes the Makefiles for the particular | |
| 106 | -system, | |
| 107 | -the <span style="font-family: Courier New,Courier,monospace;">make</span> | |
| 108 | -command compiles the source files and builds the executables, the the <span style="font-family: Courier New,Courier,monospace;">make | |
| 109 | -check</span> command performs a system test, and | |
| 110 | -the <span style="font-family: Courier New,Courier,monospace;">make | |
| 111 | -install</span> command installs the executables in the <span style="font-family: Courier New,Courier,monospace;">bin</span> | |
| 112 | -directory. | |
| 113 | -Type <span style="font-family: Courier New,Courier,monospace;">./configure</span> | |
| 114 | -and not simply <span style="font-family: Courier New,Courier,monospace;">configure</span> | |
| 115 | -to ensure that the | |
| 116 | -configuration | |
| 117 | -script in the current directory is run and not some other system-wide | |
| 118 | -configuration | |
| 119 | -script. </p> | |
| 120 | -<br> | |
| 121 | -<p style="padding-top: 0pt; text-align: left;" class="Body">By | |
| 122 | -default, the install directory is set to | |
| 123 | -<span style="font-family: Courier New,Courier,monospace;">/usr/local/gamma</span>. | |
| 124 | -To change the install directory, an optional <span style="font-family: Courier New,Courier,monospace;">--prefix</span> | |
| 125 | -argument should be | |
| 126 | -given, | |
| 127 | -for example: | |
| 128 | -</p> | |
| 129 | -<p class="Code"> | |
| 130 | -> ./configure --prefix=/usr/local | |
| 131 | -</p> | |
| 132 | -<p style="padding-top: 0pt; text-align: left;" class="Body">A | |
| 133 | -full list of configuration options can be found using</p> | |
| 134 | -<p class="Code"> | |
| 135 | -> ./configure --help | |
| 136 | -</p> | |
| 137 | -<br> | |
| 138 | -<br> | |
| 139 | -<p class="Title"><a name="Step_4:_Setup_ctatools_environment"></a>Step | |
| 140 | -4: Setup ctools environment</p> | |
| 141 | -<p style="padding-top: 0pt; text-align: left;" class="Body">If ctools have been installed in the same directory as GammaLib (which is | |
| 142 | -the default configuration for ctools) you are done.</p> | |
| 143 | -<br> | |
| 144 | -<p style="padding-top: 0pt; text-align: left;" class="Body">If | |
| 145 | -you decided to install ctools in a different directory, you have to | |
| 146 | -configure ctools by | |
| 147 | -setting up some | |
| 148 | -environment variables. | |
| 149 | -This will be done automatically by an initialisation script that is | |
| 150 | -found in the <span style="font-family: Courier New,Courier,monospace;">bin</span> | |
| 151 | -directory of the directory into which ctools has been installed. | |
| 152 | -Assuming that you have installed ctools into | |
| 153 | -<span style="font-family: Courier New,Courier,monospace;">/usr/local/ctools</span> | |
| 154 | -you need to add the following to your <span style="font-family: Courier New,Courier,monospace;">$HOME/.bashrc</span> | |
| 155 | -or <span style="font-family: Courier New,Courier,monospace;">$HOME/.profile</span> | |
| 156 | -script on a Linux | |
| 157 | -machine: | |
| 158 | -</p> | |
| 159 | -<p class="Code">export CTOOLS=/usr/local/ctools<br> | |
| 160 | -source $CTOOLS/bin/ctools-init.sh | |
| 161 | -</p> | |
| 162 | -<p style="padding-top: 0pt; text-align: left;" class="Body">If | |
| 163 | -you use C shell or a variant then add the following to your | |
| 164 | -<span style="font-family: Courier New,Courier,monospace;">$HOME/.cshrc</span> | |
| 165 | -or <span style="font-family: Courier New,Courier,monospace;">$HOME/.tcshrc</span> | |
| 166 | -script: | |
| 167 | -</p> | |
| 168 | -<p class="Code">setenv CTOOLS /usr/local/ctools<br> | |
| 169 | -source $CTOOLS/bin/ctools-init.csh </p> | |
| 170 | -<br> | |
| 171 | -<br> | |
| 172 | -<p class="Title"><a name="Known_Problems"></a>Known | |
| 173 | -Problems</p> | |
| 174 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 175 | -<span style="font-weight: bold; font-size: 14px;">Python | |
| 176 | -support</span> | |
| 177 | -<br> | |
| 178 | -ctools comes with Python wrappers so that all classes can be | |
| 179 | -directly used from Python. To compile-in Python support, ctools | |
| 180 | -needs the <span style="font-family: Courier New,Courier,monospace;">Python.h</span> | |
| 181 | -header file, which on many distributions is not | |
| 182 | -installed by default. | |
| 183 | -To make <span style="font-family: Courier New,Courier,monospace;">Python.h</span> | |
| 184 | -available, install the Python developer package | |
| 185 | -in your distribution. Otherwise you will not be able to use ctools | |
| 186 | -from Python. | |
| 187 | -</p> | |
| 188 | -<br> | |
| 189 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 190 | -<span style="font-weight: bold; font-size: 14px;">Mac OS X</span> | |
| 191 | -<br> | |
| 192 | -The Python development package is not installed as default on Mac | |
| 193 | -OS X, and consequently, the <span style="font-family: Courier New,Courier,monospace;">Python.h</span> | |
| 194 | -header file is missing that is | |
| 195 | -needed to compile in the Python bindings. The <span style="font-family: Courier New,Courier,monospace;">configure</span> | |
| 196 | -script recognises this fact and adjust the build procedure accordingly, | |
| 197 | -but you | |
| 198 | -will not be able to use ctools from Python. So better install the | |
| 199 | -Python development package before installing ctools.</p> | |
| 200 | -<p style="padding-top: 0pt; text-align: left;" class="Body">It | |
| 201 | -was also reported that adding of <span style="font-family: Courier New,Courier,monospace;">export | |
| 202 | -MACOSX_DEPLOYMENT_TARGET=10.6</span> to the <span style="font-family: Courier New,Courier,monospace;">.bashrc</span> | |
| 203 | -file was necessary one | |
| 204 | -some Mac OS X 10.6 installations to make ctools compile.<br> | |
| 205 | -</p> | |
| 206 | -<br> | |
| 207 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><span style="font-weight: bold; font-size: 14px;">Solaris</span> | |
| 208 | -<br> | |
| 209 | -Although ctools builds on Solaris using the Sun compiler, there are | |
| 210 | -problems with global symbols in shared libraries and exception | |
| 211 | -catching, which prevents the FITS interface to work | |
| 212 | -correctly. ctools has however been built and tested successfully | |
| 213 | -using the GNU | |
| 214 | -compiler, and this is the only build method that is currently | |
| 215 | -supported. | |
| 216 | -Problems have also been encountered when compiling cfitsio versions | |
| 217 | -more recent than 3.250. The problems have been reported to the cfitsio | |
| 218 | -developer team, and are likely to be solved in the future. For the | |
| 219 | -time being, it is recommended to use cfitsio version 3.250 on Solaris. | |
| 220 | -</p> | |
| 221 | -<br> | |
| 222 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><span style="font-weight: bold; font-size: 14px;">OpenSolaris</span> | |
| 223 | -<br> | |
| 224 | -On OpenSolaris, the same problems concerning the SunStudio compiler | |
| 225 | -occur as for Solaris, and also here, the GNU compiler is the | |
| 226 | -recommended | |
| 227 | -tool to build ctools. Also here, cfitsio version 3.250 is the | |
| 228 | -recommended library as more recent version feature relocation | |
| 229 | -problems. ctools has been tested using gcc 4.3.2 on OpenSolaris | |
| 230 | -2009.06. Make sure to create the symbolic links | |
| 231 | -</p> | |
| 232 | -<p class="Code">ln -s /usr/bin/gcc4.3.2 /usr/bin/gcc<br> | |
| 233 | -ln -s /usr/bin/g++4.3.2 /usr/bin/g++</p> | |
| 234 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 235 | -which are not there by default to avoid excess warnings during | |
| 236 | -compilation. | |
| 237 | -</p> | |
| 238 | -<br> | |
| 239 | -<br> | |
| 240 | -<p class="Title"><a name="Getting_Help"></a>Getting | |
| 241 | -Help</p> | |
| 242 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 243 | -Any questions, bug reports, or suggested enhancements related to ctools should be submitted to | |
| 244 | -<a href="mailto:knodlseder@cesr.fr">knodlseder@cesr.fr</a>. | |
| 245 | -</p> | |
| 246 | -<br> | |
| 247 | -<hr size="1"> | |
| 248 | -<p class="Body">Last updated by: Jürgen | |
| 249 | -Knödlseder (4 April 2012) | |
| 250 | -</p> | |
| 251 | -</body></html> | |
| 252 | 0 | \ No newline at end of file |
doc/html/quickstart.html deleted
| 1 | -<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"> | |
| 2 | -<html><head> | |
| 3 | -<meta http-equiv="content-type" content="text/html; charset=ISO-8859-1"><title>Getting Started</title> | |
| 4 | - | |
| 5 | -<link href="ctools.css" media="screen,print" type="text/css" rel="stylesheet"> | |
| 6 | -</head> | |
| 7 | -<body> | |
| 8 | -<p class="Header">Getting started</p> | |
| 9 | -<p class="Title">First steps with ctools</p> | |
| 10 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 11 | -This page gives provides a step-by-step description of the CTA data | |
| 12 | -analysis procedure implemented in ctools. In this example we will | |
| 13 | -simulate data for the Crab nebula, observed on-axis, as a single source | |
| 14 | -in the CTA field of view.</p> | |
| 15 | -<br> | |
| 16 | -<p style="padding-top: 0pt; text-align: left;" class="Body">The | |
| 17 | -example consists of the following analysis steps:</p> | |
| 18 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><a href="#Simulating_CTA_data">Simulating CTA data</a></p> | |
| 19 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><a href="#Binning_CTA_data">Binning CTA data</a></p> | |
| 20 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><a href="#Fitting_CTA_data">Fitting CTA data</a></p> | |
| 21 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><a href="#Doing_an_unbinned_analysis">Doing an unbinned analysis</a></p> | |
| 22 | -<br> | |
| 23 | -<p style="padding-top: 0pt; text-align: left;" class="Body">And | |
| 24 | -please read also: | |
| 25 | -</p> | |
| 26 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><a href="#Some_final_remarks">Some final remarks</a></p> | |
| 27 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><a href="#Getting_Help">Getting Help</a></p> | |
| 28 | -<br> | |
| 29 | -<br> | |
| 30 | -<p class="Title"><a name="Simulating_CTA_data"></a>Simulating | |
| 31 | -CTA data</p> | |
| 32 | -<p style="padding-top: 0pt; text-align: left;" class="Body">CTA | |
| 33 | -data are simulated by the executable <span style="font-family: Courier New,Courier,monospace;">ctobssim</span>. | |
| 34 | -To start the executable, type <span style="font-family: Courier New,Courier,monospace;">ctobssim</span> | |
| 35 | -at the console prompt (which is denoted by <span style="font-family: Courier New,Courier,monospace;">></span>). | |
| 36 | -It will query for a number of parameters. The example shown below will | |
| 37 | -work on a system that has been installed in the default directory <span style="font-family: Courier New,Courier,monospace;">/usr/local/gamma</span>. | |
| 38 | -If ctools have been installed in another directory, please adjust the | |
| 39 | -parameter names accordingly. | |
| 40 | -</p> | |
| 41 | -<p class="Code"> | |
| 42 | -> ctobssim<br> | |
| 43 | -Model [test/data/crab_alone.xml] /usr/local/gamma/share/models/crab.xml | |
| 44 | -<br> | |
| 45 | -Output event data file [events.fits] <br> | |
| 46 | -Calibration database [test/irf] /usr/local/gamma/share/caldb/cta<br> | |
| 47 | -Instrument response function [kb_E_50h_v3] <br> | |
| 48 | -RA for ROI centre (degrees) (0-360) [83.63] <br> | |
| 49 | -Dec for ROI centre (degrees) (-90-90) [22.01] <br> | |
| 50 | -Radius of ROI (degrees) (0-180) [10.0] <br> | |
| 51 | -Start time (MET in s) (0) [0.0] <br> | |
| 52 | -End time (MET in s) (0) [1800.0] <br> | |
| 53 | -Lower energy limit (TeV) (0) [0.1] <br> | |
| 54 | -Upper energy limit (TeV) (0) [100.0] <br> | |
| 55 | -</p> | |
| 56 | -<p style="padding-top: 0pt; text-align: left;" class="Body">Each | |
| 57 | -line represents a query for one parameter. The line starts with a short | |
| 58 | -description of the parameter, followed by the default parameter value | |
| 59 | -proposed by <span style="font-family: Courier New,Courier,monospace;">ctobssim</span> | |
| 60 | -in squared brackets <span style="font-family: Courier New,Courier,monospace;">[ ]</span>. | |
| 61 | -If no parameter is entered (which is the case for the majority of | |
| 62 | -parameters shown here), the default parameter will be used.</p> | |
| 63 | -<br> | |
| 64 | -<p style="padding-top: 0pt; text-align: left;" class="Body">The | |
| 65 | -CTA instrument properties (effective area, PSF width) are taken for the | |
| 66 | -moment from Konrad Bernlöhrs sensitivity files which are | |
| 67 | -stored in <span style="font-family: Courier New,Courier,monospace;">/usr/local/gamma/share/caldb/cta</span>. | |
| 68 | -To test another CTA configuration you may adapt the instrument response | |
| 69 | -function name accordingly (see the <a href="#Some_final_remarks">final | |
| 70 | -remarks</a>).</p> | |
| 71 | -<br> | |
| 72 | -<p style="padding-top: 0pt; text-align: left;" class="Body">Events | |
| 73 | -are simulated based on the instrument properties and based on a source | |
| 74 | -and background model. Only events that fall within the | |
| 75 | -specified region of interest (ROI), defined as a circle around a sky | |
| 76 | -position in Right Ascension and Declination (in degrees), will be | |
| 77 | -stored in the output event data file. The duration of the simulation is | |
| 78 | -taken here to 30 minutes (or 1800 seconds). Events are simulated for | |
| 79 | -energies between 0.1 and 100 TeV.</p> | |
| 80 | -<br> | |
| 81 | -<p style="padding-top: 0pt; text-align: left;" class="Body">The | |
| 82 | -source and background model is defined by the XML file <span style="font-family: Courier New,Courier,monospace;">/usr/local/gamma/share/models/crab.xml</span>: | |
| 83 | -</p> | |
| 84 | -<p class="Code"> | |
| 85 | -<?xml version="1.0" standalone="no"?><br> | |
| 86 | -<source_library title="source library"><br> | |
| 87 | - <source name="MAGIC Crab (Albert et al. 2008, ApJ, | |
| 88 | -674, 1037)" type="PointSource"><br> | |
| 89 | - <spectrum type="PowerLaw"><br> | |
| 90 | - | |
| 91 | -<parameter name="Prefactor" scale="1e-16" value="5.7" | |
| 92 | -min="1e-07" max="1000.0" free="1"/><br> | |
| 93 | - | |
| 94 | -<parameter name="Index" | |
| 95 | -scale="-1" value="2.48" | |
| 96 | -min="0.0" max="+5.0" | |
| 97 | -free="1"/><br> | |
| 98 | - | |
| 99 | -<parameter name="Scale" | |
| 100 | -scale="1e6" value="0.3" | |
| 101 | -min="0.01" max="1000.0" free="0"/><br> | |
| 102 | - </spectrum><br> | |
| 103 | - <spatialModel | |
| 104 | -type="SkyDirFunction"><br> | |
| 105 | - <parameter | |
| 106 | -name="RA" scale="1.0" value="83.6331" min="-360" max="360" | |
| 107 | -free="0"/><br> | |
| 108 | - <parameter | |
| 109 | -name="DEC" scale="1.0" value="22.0145" min="-90" | |
| 110 | -max="90" free="0"/><br> | |
| 111 | - </spatialModel><br> | |
| 112 | - </source><br> | |
| 113 | - <source name="Background" type="RadialAcceptance" | |
| 114 | -instrument="CTA"><br> | |
| 115 | - <spectrum type="PowerLaw"><br> | |
| 116 | - | |
| 117 | -<parameter name="Prefactor" scale="1e-6" value="61.8" | |
| 118 | -min="0.0" max="1000.0" free="1"/><br> | |
| 119 | - | |
| 120 | -<parameter name="Index" | |
| 121 | -scale="-1" value="1.85" min="0.0" | |
| 122 | -max="+5.0" free="1"/><br> | |
| 123 | - | |
| 124 | -<parameter name="Scale" | |
| 125 | -scale="1e6" value="1.0" min="0.01" max="1000.0" | |
| 126 | -free="0"/><br> | |
| 127 | - </spectrum><br> | |
| 128 | - <radialModel | |
| 129 | -type="Gaussian"><br> | |
| 130 | - | |
| 131 | -<parameter name="Sigma" scale="1.0" value="3.0" min="0.01" | |
| 132 | -max="10.0" free="1"/><br> | |
| 133 | - </radialModel><br> | |
| 134 | - </source><br> | |
| 135 | -</source_library> | |
| 136 | -</p> | |
| 137 | -<p style="padding-top: 0pt; text-align: left;" class="Body">The | |
| 138 | -model consists of a source library that contains 2 | |
| 139 | -"sources": the Crab nebula and an instrumental background | |
| 140 | -model.</p> | |
| 141 | -<br> | |
| 142 | -<p style="padding-top: 0pt; text-align: left;" class="Body">The | |
| 143 | -Crab nebula is modelled by a factorized sky model that has a spectral | |
| 144 | -and a spatial component (tags <span style="font-family: Courier New,Courier,monospace;"><spectrum></span> | |
| 145 | -and <span style="font-family: Courier New,Courier,monospace;"><spatialModel></span>, | |
| 146 | -respectively). The spectrum is modelled by a power law, which is | |
| 147 | -defined by 3 parameters: the <span style="font-family: Courier New,Courier,monospace;">Prefactor</span>, | |
| 148 | -the <span style="font-family: Courier New,Courier,monospace;">Index</span> | |
| 149 | -and the <span style="font-family: Courier New,Courier,monospace;">Scale</span>. | |
| 150 | -The spatial model has 2 parameters: Right Ascension in degrees (<span style="font-family: Courier New,Courier,monospace;">RA</span>), | |
| 151 | -and Declination in degrees (<span style="font-family: Courier New,Courier,monospace;">DEC</span>). | |
| 152 | -Each parameter has a <span style="font-family: Courier New,Courier,monospace;">value</span> | |
| 153 | -and a <span style="font-family: Courier New,Courier,monospace;">scale</span> | |
| 154 | -factor, the real value of the parameter being the product <span style="font-family: Courier New,Courier,monospace;">value * | |
| 155 | -scale</span>. Typically, <span style="font-family: Courier New,Courier,monospace;">scale</span> | |
| 156 | -is chosen so that <span style="font-family: Courier New,Courier,monospace;">value</span> | |
| 157 | -is of the order of 1 (this is relevant for model fitting later). In | |
| 158 | -addition, <span style="font-family: Courier New,Courier,monospace;">value</span> | |
| 159 | -is bound by a minimum (<span style="font-family: Courier New,Courier,monospace;">min</span>) | |
| 160 | -and maximum (<span style="font-family: Courier New,Courier,monospace;">max</span>) | |
| 161 | -value, and a parameter may be free (<span style="font-family: Courier New,Courier,monospace;">free="1"</span>) | |
| 162 | -or fixed (<span style="font-family: Courier New,Courier,monospace;">free="0"</span>). | |
| 163 | -The <span style="font-family: Courier New,Courier,monospace;">min</span>, | |
| 164 | -<span style="font-family: Courier New,Courier,monospace;">max</span>, | |
| 165 | -and <span style="font-family: Courier New,Courier,monospace;">free</span> | |
| 166 | -attributes are not relevant here for the simulations, but they will be | |
| 167 | -important for the model fitting later.</p> | |
| 168 | -<br> | |
| 169 | -<p style="padding-top: 0pt; text-align: left;" class="Body">The | |
| 170 | -spectral intensity <span style="font-family: Courier New,Courier,monospace;">I(E)</span> | |
| 171 | -(in units of photons/cm2/s/MeV) of the power law is given by </p> | |
| 172 | -<p style="padding-top: 0pt; text-align: center; font-family: Courier New,Courier,monospace;" class="Body">I(E) = Prefactor * (E/Scale)^Index</p> | |
| 173 | -<p style="padding-top: 0pt; text-align: left;" class="Body">where | |
| 174 | -the <span style="font-family: Courier New,Courier,monospace;">E</span> | |
| 175 | -and <span style="font-family: Courier New,Courier,monospace;">Scale</span> | |
| 176 | -are given in MeV.</p> | |
| 177 | -<br> | |
| 178 | -<p style="padding-top: 0pt; text-align: left;" class="Body">The | |
| 179 | -instrumental background of CTA is modelled by a factorized data model | |
| 180 | -that has a spectral and a radial component (tags <span style="font-family: Courier New,Courier,monospace;"><spectrum></span> | |
| 181 | -and <span style="font-family: Courier New,Courier,monospace;"><radialModel></span>, | |
| 182 | -respectively). The spectral component describes the <span style="font-weight: bold;">on-axis background | |
| 183 | -counting rate</span> of CTA as function of energy in units of | |
| 184 | -counts/s/sr/TeV. | |
| 185 | -The radial component describes the variation of the background rate | |
| 186 | -with <span style="font-weight: bold;">offset angle squared</span>, | |
| 187 | -(i.e. square of the offset angle with respect to the pointing | |
| 188 | -direction) which is modelled here by a Gaussian. The only parameter of | |
| 189 | -the radial component is the width of the Gaussian <span style="font-family: Courier New,Courier,monospace;">Sigma</span>, | |
| 190 | -which is here set to 3 degrees squared.</p> | |
| 191 | -<br> | |
| 192 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><span style="font-family: Courier New,Courier,monospace;">ctobssim</span> | |
| 193 | -has a couple of hidden parameters, the most important one being | |
| 194 | -certainly <span style="font-family: Courier New,Courier,monospace;">seed</span>. | |
| 195 | -<span style="font-family: Courier New,Courier,monospace;">seed</span> | |
| 196 | -is an integer that specifies the seed value for the random number | |
| 197 | -generator, and changing this parameter will allow to generate | |
| 198 | -statistically independent Monte Carlo samples of CTA event data. To | |
| 199 | -use for example a seed value of 41 one should type: | |
| 200 | -</p> | |
| 201 | -<p class="Code">> ctobssim seed=41<br> | |
| 202 | -</p> | |
| 203 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 204 | -<span style="font-family: Courier New,Courier,monospace;">ctobssim</span> | |
| 205 | -will write 2 files in the working directory: <span style="font-family: Courier New,Courier,monospace;">events.fits</span> | |
| 206 | -and <span style="font-family: Courier New,Courier,monospace;">ctobssim.log</span>. | |
| 207 | -The first file contains the simulated events in FITS format and can be | |
| 208 | -inspected using <a href="http://heasarc.gsfc.nasa.gov/docs/software/ftools/fv/" target="_blank"><span style="font-family: Courier New,Courier,monospace;">fv</span></a> | |
| 209 | -or <a href="http://hea-www.harvard.edu/RD/ds9/" target="_blank"><span style="font-family: Courier New,Courier,monospace;">ds9</span></a>. | |
| 210 | -The FITS file will contain 3 extensions: an empty primary image, a | |
| 211 | -binary table named <span style="font-family: Courier New,Courier,monospace;">EVENTS</span> | |
| 212 | -that holds the events (one row per event), and a binary table named <span style="font-family: Courier New,Courier,monospace;">GTI</span> | |
| 213 | -holding the Good Time Interval (for the moment a single row with 2 | |
| 214 | -columns providing the start and the stop time of the simulated time | |
| 215 | -interval).</p> | |
| 216 | -<br> | |
| 217 | -<p style="padding-top: 0pt; text-align: left;" class="Body">The | |
| 218 | -second file produced by <span style="font-family: Courier New,Courier,monospace;">ctobssim</span> | |
| 219 | -is a human readable log file that contains information about | |
| 220 | -the job execution. As example, the last lines from this file are shown | |
| 221 | -here:</p> | |
| 222 | -<p class="Code">2011-01-31T23:50:39: +==================+<br> | |
| 223 | -2011-01-31T23:50:39: | Simulate photons |<br> | |
| 224 | -2011-01-31T23:50:39: +==================+<br> | |
| 225 | -2011-01-31T23:50:39: Simulated photons .........: 33943<br> | |
| 226 | -2011-01-31T23:50:39: Simulation area ...........: 3.2e+10 cm2<br> | |
| 227 | -2011-01-31T23:50:39: Simulation cone ...........: RA=83.63 | |
| 228 | -deg, Dec=22.01deg, r=10 deg<br> | |
| 229 | -2011-01-31T23:50:39: Time interval .............: 0 s - 1800 s<br> | |
| 230 | -2011-01-31T23:50:39: Energy range ..............: 100 GeV - | |
| 231 | -100 TeV<br> | |
| 232 | -2011-01-31T23:50:39: Random number seed ........: 1<br> | |
| 233 | -2011-01-31T23:50:39: <br> | |
| 234 | -2011-01-31T23:50:39: +=================+<br> | |
| 235 | -2011-01-31T23:50:39: | Simulate events |<br> | |
| 236 | -2011-01-31T23:50:39: +=================+<br> | |
| 237 | -2011-01-31T23:50:39: MC source events ..........: 649<br> | |
| 238 | -2011-01-31T23:50:39: MC background events ......: 3407<br> | |
| 239 | -2011-01-31T23:50:39: <br> | |
| 240 | -2011-01-31T23:50:39: +==========================+<br> | |
| 241 | -2011-01-31T23:50:39: | Save events in FITS file |<br> | |
| 242 | -2011-01-31T23:50:39: +==========================+<br> | |
| 243 | -2011-01-31T23:50:39: <br> | |
| 244 | -2011-01-31T23:50:39: Application "ctobssim" terminated after consuming | |
| 245 | -14 seconds of CPU time.<br> | |
| 246 | -<br> | |
| 247 | -</p> | |
| 248 | -<p style="padding-top: 0pt; text-align: left;" class="Body">Each | |
| 249 | -line starts with the UTC time at which the line has been written. In | |
| 250 | -this | |
| 251 | -run, 33943 Crab photons have been thrown over an area of 3.2 square | |
| 252 | -kilometres during a time interval of 1800 seconds. 649 of these photons | |
| 253 | -have been registered by CTA as events. In the same time | |
| 254 | -interval, 3407 background events have been registred by CTA.</p> | |
| 255 | -<br> | |
| 256 | -<br> | |
| 257 | -<p class="Title"><a name="Binning_CTA_data"></a>Binning | |
| 258 | -CTA data</p> | |
| 259 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 260 | -As next analysis step the data are binned in a counts map using the | |
| 261 | -executable <span style="font-family: Courier New,Courier,monospace;">ctbin</span>. | |
| 262 | -A counts map is a 3 dimensional data cube, spanned by Right Ascension | |
| 263 | -(or Galactic longitude), Declination (or Galactic latitude), and the | |
| 264 | -logarithm (base 10) of energy.</p> | |
| 265 | -<br> | |
| 266 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><span style="font-family: Courier New,Courier,monospace;">ctbin</span> | |
| 267 | -is executed by typing:</p> | |
| 268 | -<p class="Code"> | |
| 269 | -> ctbin<br> | |
| 270 | -Event data file name [test/data/crab_events.fits.gz] events.fits<br> | |
| 271 | -Output file name [cntmap.fits] <br> | |
| 272 | -Start value for first energy bin in TeV [0.1] <br> | |
| 273 | -Stop value for last energy bin in TeV [100.0] <br> | |
| 274 | -Number of energy bins [20] <br> | |
| 275 | -Projection method e.g. AIT|ARC|CAR|GLS|MER|NCP|SIN|STG|TAN [CAR] <br> | |
| 276 | -Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] <br> | |
| 277 | -First coordinate of image center in degrees (RA or galactic l) [83.63] <br> | |
| 278 | -Second coordinate of image center in degrees (DEC or galactic b) | |
| 279 | -[22.01] <br> | |
| 280 | -Image scale (in degrees/pixel) [0.02] <br> | |
| 281 | -Size of the X axis in pixels [200] <br> | |
| 282 | -Size of the Y axis in pixels [200] </p> | |
| 283 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 284 | -In this example we adjust the event data file name and accept all the | |
| 285 | -remaining parameter defaults as they perfectly | |
| 286 | -satisfy our needs. The counts map will be centred on the location of | |
| 287 | -the Crab (Right Ascension 83.63 degrees, Declination 22.01 degrees) and | |
| 288 | -will be aligned in celestial coordinates. A cartesian projection has | |
| 289 | -been selected (<span style="font-weight: bold;">note that | |
| 290 | -the cartesian projection is the only valid project available for the | |
| 291 | -moment</span>). The counts map has 200 x 200 | |
| 292 | -spatial pixels of 0.02 x 0.02 degrees in size, hence it covers a total | |
| 293 | -area of 4 x 4 degrees.</p> | |
| 294 | -<br> | |
| 295 | -<p style="padding-top: 0pt; text-align: left;" class="Body">A | |
| 296 | -set of 20 counts maps will be created, which are logarithmically spaced | |
| 297 | -in energy, and which cover the energy range from 0.1 TeV to 100 TeV. In | |
| 298 | -this example, the counts map will be saved as <span style="font-family: Courier New,Courier,monospace;">cntmap.fits</span> | |
| 299 | -in the | |
| 300 | -working directory. | |
| 301 | -In addition to the counts map, that is stored as the primary image | |
| 302 | -extension, the FITS | |
| 303 | -file also contains an extension named <span style="font-family: Courier New,Courier,monospace;">EBOUNDS</span> | |
| 304 | -that defines the energy boundaries that were used, and an extension <span style="font-family: Courier New,Courier,monospace;">GTI</span> | |
| 305 | -that defines the Good Time Intervals that have been used. The following | |
| 306 | -image shows the resulting FITS file. The <span style="font-family: Courier New,Courier,monospace;">EBOUNDS</span> | |
| 307 | -table has 20 rows, one for each energy bin, while the <span style="font-family: Courier New,Courier,monospace;">GTI</span> | |
| 308 | -table has jusr a single row, indicating the start and stop time of the | |
| 309 | -simulated data.</p> | |
| 310 | -<br> | |
| 311 | -<div style="text-align: center;"><img alt="CTA counts map FITS file" src="images/cntmap-fits.jpg"><br> | |
| 312 | -</div> | |
| 313 | -<br> | |
| 314 | -<p style="padding-top: 0pt; text-align: left;" class="Body">An | |
| 315 | -image | |
| 316 | -of the first bin, covering the energy range 100 - 141 GeV, is | |
| 317 | -shown below: | |
| 318 | -</p> | |
| 319 | -<div style="text-align: center;"><img style="width: 496px; height: 518px;" alt="CTA counts map" src="images/cntmap-map.jpg"><br> | |
| 320 | -</div> | |
| 321 | -<br> | |
| 322 | -<p style="padding-top: 0pt; text-align: left;" class="Body">For | |
| 323 | -illustrating, the last few lines of the log file <span style="font-family: Courier New,Courier,monospace;">ctbin.log</span> | |
| 324 | -are shown below: | |
| 325 | -</p> | |
| 326 | -<p class="Code">2011-02-01T09:28:05: +=========+<br> | |
| 327 | -2011-02-01T09:28:05: | Binning |<br> | |
| 328 | -2011-02-01T09:28:05: +=========+<br> | |
| 329 | -2011-02-01T09:28:05: Events in list ............: 4056<br> | |
| 330 | -2011-02-01T09:28:05: Events in map .............: 3683<br> | |
| 331 | -2011-02-01T09:28:05: Events outside map area ...: 373<br> | |
| 332 | -2011-02-01T09:28:05: Events outside energy bins : 0<br> | |
| 333 | -2011-02-01T09:28:05: <br> | |
| 334 | -2011-02-01T09:28:05: +============+<br> | |
| 335 | -2011-02-01T09:28:05: | Counts map |<br> | |
| 336 | -2011-02-01T09:28:05: +============+<br> | |
| 337 | -2011-02-01T09:28:05: === GSkymap ===<br> | |
| 338 | -2011-02-01T09:28:05: Number of pixels ..........: 40000<br> | |
| 339 | -2011-02-01T09:28:05: Number of maps ............: 20<br> | |
| 340 | -2011-02-01T09:28:05: X axis dimension ..........: 200<br> | |
| 341 | -2011-02-01T09:28:05: Y axis dimension ..........: 200<br> | |
| 342 | -2011-02-01T09:28:05: === GWcsCAR ===<br> | |
| 343 | -2011-02-01T09:28:05: Coordinate system .........: EQU<br> | |
| 344 | -2011-02-01T09:28:05: Reference coordinate ......: (83.63, | |
| 345 | -22.01) deg<br> | |
| 346 | -2011-02-01T09:28:05: Reference pixel ...........: (100.5, | |
| 347 | -100.5)<br> | |
| 348 | -2011-02-01T09:28:05: Increment at reference ....: (0.02, | |
| 349 | -0.02) deg<br> | |
| 350 | -2011-02-01T09:28:05: Coordinate of North Pole ..: (180, 0) deg<br> | |
| 351 | -2011-02-01T09:28:05: CD matrix .................: [0.02 0][0 | |
| 352 | -0.02]<br> | |
| 353 | -2011-02-01T09:28:05: Inverse CD matrix .........: [50 -0][-0 | |
| 354 | -50]<br> | |
| 355 | -2011-02-01T09:28:05: Coordinate of native pole .: (83.63, | |
| 356 | -112.01)<br> | |
| 357 | -2011-02-01T09:28:05: <br> | |
| 358 | -2011-02-01T09:28:05: Application "ctbin" terminated after consuming 10 | |
| 359 | -seconds of CPU time. | |
| 360 | -</p> | |
| 361 | -<p style="padding-top: 0pt; text-align: left;" class="Body">From | |
| 362 | -the 4056 events that have been simulated and stored in the <span style="font-family: Courier New,Courier,monospace;">events.fits</span> | |
| 363 | -file, 3683 lie within the map boundaries and are thus put into the | |
| 364 | -resulting counts map. The counts map is stored in a cartesian | |
| 365 | -projection in a <a href="http://fits.gsfc.nasa.gov/fits_wcs.html" target="_blank">World Coordinate System (WCS)</a> | |
| 366 | -compliant format. | |
| 367 | -</p> | |
| 368 | -<br> | |
| 369 | -<br> | |
| 370 | -<p class="Title"><a name="Fitting_CTA_data"></a>Fitting | |
| 371 | -CTA data</p> | |
| 372 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 373 | -Now we are ready to fit the simulated data with a model. For simplicity | |
| 374 | -we use in this example the same model that we used to simulate the data | |
| 375 | -with <span style="font-family: Courier New,Courier,monospace;">ctobssim</span>. | |
| 376 | -Model fitting is done using the executable <span style="font-family: Courier New,Courier,monospace;">ctlike</span>, | |
| 377 | -and we do the fit by typing: | |
| 378 | -</p> | |
| 379 | -<p class="Code"> | |
| 380 | -> ctlike<br> | |
| 381 | -Event list, counts map or observation definition | |
| 382 | -file [test/data/crab_cntmap.fits.gz] cntmap.fits<br> | |
| 383 | -Calibration | |
| 384 | -database [test/irf] /usr/local/gamma/share/caldb/cta<br> | |
| 385 | -Instrument response function [kb_E_50h_v3] <br> | |
| 386 | -Source model [test/data/crab.xml] /usr/local/gamma/share/models/crab.xml<br> | |
| 387 | -Source model output file [crab_results.xml] | |
| 388 | -</p> | |
| 389 | -<p style="padding-top: 0pt; text-align: left;" class="Body">Fitting | |
| 390 | -of the data is done in <span style="font-family: Courier New,Courier,monospace;">BINNED</span> | |
| 391 | -mode, which means that the events have been binned into a counts map | |
| 392 | -and the fit computes the log-likelihood function by summing over all | |
| 393 | -200 x 200 x 20 bins of the counts map. There is an alternative method, | |
| 394 | -the so called <span style="font-family: Courier New,Courier,monospace;">UNBINNED</span> | |
| 395 | -mode, where the events are not binned into a counts map and the | |
| 396 | -log-likelihood is computed directly by summing over all events. We will | |
| 397 | -explore | |
| 398 | -the <span style="font-family: Courier New,Courier,monospace;">UNBINNED</span> | |
| 399 | -mode <a href="#Doing_an_unbinned_analysis">later</a>.</p> | |
| 400 | -<br> | |
| 401 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 402 | -One of the parameters given to <span style="font-family: Courier New,Courier,monospace;">ctlike</span> | |
| 403 | -is a source model output file (we specified <span style="font-family: Courier New,Courier,monospace;">crab_results.xml</span> | |
| 404 | -in the example), and this file will be a copy of the model XML file | |
| 405 | -where the parameter values have been replaced by the fit results. In | |
| 406 | -addition, the statistical uncertainties are added for each fitted | |
| 407 | -parameter using the attribute <span style="font-family: Courier New,Courier,monospace;">error</span>. | |
| 408 | -Below we show the XML result file that has been produced by the run: | |
| 409 | -</p> | |
| 410 | -<p class="Code"><?xml version="1.0" | |
| 411 | -encoding="UTF-8" standalone="no"?><br> | |
| 412 | -<source_library title="source library"><br> | |
| 413 | - <source name="MAGIC Crab (Albert et al. 2008, ApJ, | |
| 414 | -674, 1037)" type="DiffuseSource"><br> | |
| 415 | - <spectrum type="PowerLaw"><br> | |
| 416 | - <parameter | |
| 417 | -name="Prefactor" value="5.7765" error="0.243328" scale="1e-16" | |
| 418 | -min="1e-07" max="1000" free="1" /><br> | |
| 419 | - <parameter | |
| 420 | -name="Index" value="2.49937" error="0.0330191" scale="-1" min="0" | |
| 421 | -max="5" free="1" /><br> | |
| 422 | - <parameter | |
| 423 | -name="Scale" value="0.3" scale="1e+06" min="0.01" max="1000" free="0" | |
| 424 | -/><br> | |
| 425 | - </spectrum><br> | |
| 426 | - <spatialModel | |
| 427 | -type="SkyDirFunction"><br> | |
| 428 | - <parameter | |
| 429 | -name="RA" value="83.6331" scale="1" min="-360" max="360" free="0" | |
| 430 | -/><br> | |
| 431 | - <parameter | |
| 432 | -name="DEC" value="22.0145" scale="1" min="-90" max="90" free="0" | |
| 433 | -/><br> | |
| 434 | - </spatialModel><br> | |
| 435 | - </source><br> | |
| 436 | - <source name="Background" type="RadialAcceptance" | |
| 437 | -instrument="CTA"><br> | |
| 438 | - <spectrum type="PowerLaw"><br> | |
| 439 | - <parameter | |
| 440 | -name="Prefactor" value="61.7696" error="1.89829" scale="1e-06" min="0" | |
| 441 | -max="1000" free="1" /><br> | |
| 442 | - <parameter | |
| 443 | -name="Index" value="1.87195" error="0.0166039" scale="-1" min="0" | |
| 444 | -max="5" free="1" /><br> | |
| 445 | - <parameter | |
| 446 | -name="Scale" value="1" scale="1e+06" min="0.01" max="1000" free="0" | |
| 447 | -/><br> | |
| 448 | - </spectrum><br> | |
| 449 | - <radialModel | |
| 450 | -type="Gaussian"><br> | |
| 451 | - <parameter | |
| 452 | -name="Sigma" value="3.00467" error="0.071261" scale="1" min="0.01" | |
| 453 | -max="10" free="1" /><br> | |
| 454 | - </radialModel><br> | |
| 455 | - </source><br> | |
| 456 | -</source_library> | |
| 457 | -</p> | |
| 458 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 459 | -In this example, the <span style="font-family: Courier New,Courier,monospace;">Prefactor</span> | |
| 460 | -and <span style="font-family: Courier New,Courier,monospace;">Index</span> | |
| 461 | -of the spectral model for the Crab as well as the <span style="font-family: Courier New,Courier,monospace;">Prefactor</span>, | |
| 462 | -<span style="font-family: Courier New,Courier,monospace;">Index</span> | |
| 463 | -and <span style="font-family: Courier New,Courier,monospace;">Sigma</span> | |
| 464 | -parameter of the radial acceptance model have been fitted (all | |
| 465 | -parameters having the attribute <span style="font-family: Courier New,Courier,monospace;">free="1"</span> | |
| 466 | -are fitted). | |
| 467 | -</p> | |
| 468 | -<br> | |
| 469 | -<p style="padding-top: 0pt; text-align: left;" class="Body">To | |
| 470 | -get more details about the model fitting one can inspect the log file. | |
| 471 | -Below | |
| 472 | -we show the last lines of the <span style="font-family: Courier New,Courier,monospace;">ctlike.log</span> | |
| 473 | -log file that has been produced by this run: | |
| 474 | -</p> | |
| 475 | -<p class="Code">2011-02-01T09:53:17: | |
| 476 | -+=================================+<br> | |
| 477 | -2011-02-01T09:53:17: | Maximum likelihood optimisation |<br> | |
| 478 | -2011-02-01T09:53:17: +=================================+<br> | |
| 479 | -2011-02-01T09:53:20: Initial iteration: func=18092.1, | |
| 480 | -Lambda=0.001<br> | |
| 481 | -2011-02-01T09:53:24: Iteration 1: func=18090, Lambda=0.0001, | |
| 482 | -delta=2.08792<br> | |
| 483 | -2011-02-01T09:53:27: Iteration 2: func=18090, Lambda=1e-05, | |
| 484 | -delta=0.0021398<br> | |
| 485 | -2011-02-01T09:53:30: Iteration 3: func=18090, Lambda=1e-06, | |
| 486 | -delta=4.80215e-06<br> | |
| 487 | -2011-02-01T09:53:34: Iteration 4: func=18090, Lambda=1e-07, | |
| 488 | -delta=1.7495e-08<br> | |
| 489 | -2011-02-01T09:53:37: <br> | |
| 490 | -2011-02-01T09:53:37: +=========================================+<br> | |
| 491 | -2011-02-01T09:53:37: | Maximum likelihood optimisation results |<br> | |
| 492 | -2011-02-01T09:53:37: +=========================================+<br> | |
| 493 | -2011-02-01T09:53:37: === GOptimizerLM ===<br> | |
| 494 | -2011-02-01T09:53:37: Optimized function value ..: 18090<br> | |
| 495 | -2011-02-01T09:53:37: Absolute precision ........: 1e-06<br> | |
| 496 | -2011-02-01T09:53:37: Optimization status .......: converged<br> | |
| 497 | -2011-02-01T09:53:37: Number of parameters ......: 11<br> | |
| 498 | -2011-02-01T09:53:37: Number of free parameters .: 5<br> | |
| 499 | -2011-02-01T09:53:37: Number of iterations ......: 4<br> | |
| 500 | -2011-02-01T09:53:37: Lambda ....................: 1e-07<br> | |
| 501 | -2011-02-01T09:53:37: Maximum log likelihood ....: -18090<br> | |
| 502 | -2011-02-01T09:53:37: Npred .....................: 3683<br> | |
| 503 | -2011-02-01T09:53:37: === GModels ===<br> | |
| 504 | -2011-02-01T09:53:37: Number of models ..........: 2<br> | |
| 505 | -2011-02-01T09:53:37: Number of parameters ......: 11<br> | |
| 506 | -2011-02-01T09:53:37: === GModelPointSource ===<br> | |
| 507 | -2011-02-01T09:53:37: Name ......................: MAGIC Crab | |
| 508 | -(Albert et al. 2008, ApJ, 674, 1037)<br> | |
| 509 | -2011-02-01T09:53:37: Instruments ...............: all<br> | |
| 510 | -2011-02-01T09:53:37: Number of parameters ......: 6<br> | |
| 511 | -2011-02-01T09:53:37: Number of spatial par's ...: 2<br> | |
| 512 | -2011-02-01T09:53:37: RA .......................: | |
| 513 | -83.6331 [-360,360] deg (fixed,scale=1)<br> | |
| 514 | -2011-02-01T09:53:37: DEC ......................: | |
| 515 | -22.0145 [-90,90] deg (fixed,scale=1)<br> | |
| 516 | -2011-02-01T09:53:37: Number of spectral par's ..: 3<br> | |
| 517 | -2011-02-01T09:53:37: Prefactor ................: | |
| 518 | -5.7765e-16 +/- 2.43328e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16)<br> | |
| 519 | -2011-02-01T09:53:37: Index ....................: | |
| 520 | --2.49937 +/- 0.0330191 [-0,-5] (free,scale=-1)<br> | |
| 521 | -2011-02-01T09:53:37: PivotEnergy ..............: | |
| 522 | -300000 [10000,1e+09] MeV (fixed,scale=1e+06)<br> | |
| 523 | -2011-02-01T09:53:37: Number of temporal par's ..: 1<br> | |
| 524 | -2011-02-01T09:53:37: Constant .................: 1 | |
| 525 | -(relative value) (fixed,scale=1)<br> | |
| 526 | -2011-02-01T09:53:37: === GCTAModelRadialAcceptance ===<br> | |
| 527 | -2011-02-01T09:53:37: Name ......................: Background<br> | |
| 528 | -2011-02-01T09:53:37: Instruments ...............: CTA<br> | |
| 529 | -2011-02-01T09:53:37: Model type ................: "Gaussian" | |
| 530 | -* "PowerLaw" * "Constant"<br> | |
| 531 | -2011-02-01T09:53:37: Number of parameters ......: 5<br> | |
| 532 | -2011-02-01T09:53:37: Number of radial par's ....: 1<br> | |
| 533 | -2011-02-01T09:53:37: Sigma ....................: | |
| 534 | -3.00467 +/- 0.071261 [0.01,10] deg2 (free,scale=1)<br> | |
| 535 | -2011-02-01T09:53:37: Number of spectral par's ..: 3<br> | |
| 536 | -2011-02-01T09:53:37: Prefactor ................: | |
| 537 | -6.17696e-05 +/- 1.89829e-06 [0,0.001] ph/cm2/s/MeV (free,scale=1e-06)<br> | |
| 538 | -2011-02-01T09:53:37: Index ....................: | |
| 539 | --1.87195 +/- 0.0166039 [-0,-5] (free,scale=-1)<br> | |
| 540 | -2011-02-01T09:53:37: PivotEnergy ..............: | |
| 541 | -1e+06 [10000,1e+09] MeV (fixed,scale=1e+06)<br> | |
| 542 | -2011-02-01T09:53:37: Number of temporal par's ..: 1<br> | |
| 543 | -2011-02-01T09:53:37: Constant .................: 1 | |
| 544 | -(relative value) (fixed,scale=1)<br> | |
| 545 | -2011-02-01T09:53:37: <br> | |
| 546 | -2011-02-01T09:53:37: Application "ctlike" terminated after consuming 41 | |
| 547 | -seconds of CPU time. | |
| 548 | -</p> | |
| 549 | -<p style="padding-top: 0pt; text-align: left;" class="Body">The | |
| 550 | -maximum likelihood optimizer required 4 iterations to converge. This is | |
| 551 | -pretty fast, but recall that we used the same model file for the | |
| 552 | -simulation and for fitting, hence the initial parameter values are | |
| 553 | -already very close to the best fitting values. To see the impact of the | |
| 554 | -initial parameters on the fit result, you may re-run <span style="font-family: Courier New,Courier,monospace;">ctlike</span> | |
| 555 | -using another copy of the model XML file where the <span style="font-family: Courier New,Courier,monospace;">value</span> | |
| 556 | -attributes of the parameters that should be fitted are changed. You | |
| 557 | -will see that the optimizer requires a couple of more iterations, but | |
| 558 | -it should converge to the same solution (provided that the initial | |
| 559 | -values are not too far of the best fitting values). | |
| 560 | -</p> | |
| 561 | -<br> | |
| 562 | -<br> | |
| 563 | -<p class="Title"><a name="Doing_an_unbinned_analysis"></a>Doing | |
| 564 | -an unbinned analysis</p> | |
| 565 | -<p style="padding-top: 0pt; text-align: left;" class="Body">As | |
| 566 | -gamma-ray events are rare, the counts maps generated by <span style="font-family: Courier New,Courier,monospace;">ctbin</span> | |
| 567 | -will in general be sparse, having many empty pixels, in | |
| 568 | -particular at high energies. An alternative analysis technique consists | |
| 569 | -of working directly on the event list without binning the events in a | |
| 570 | -counts map. We | |
| 571 | -will see the benefit of such an analysis later once we re-run <span style="font-family: Courier New,Courier,monospace;">ctlike</span> | |
| 572 | -in unbinned mode.</p> | |
| 573 | -<br> | |
| 574 | -<p style="padding-top: 0pt; text-align: left;" class="Body">For | |
| 575 | -unbinned analysis we have first to define the data | |
| 576 | -space region over which the analysis is done. This is similiar to the <span style="font-family: Courier New,Courier,monospace;">ctbin</span> | |
| 577 | -step in binned analysis where we define the size of the counts map, the | |
| 578 | -energy range, and the time interval. For unbinned analysis we have no | |
| 579 | -such thing as a counts map, but we have to define over which region of | |
| 580 | -the data space the selected events are spread (because we have to | |
| 581 | -integrate | |
| 582 | -over this region to compute the total number of predicted events in the | |
| 583 | -data space that we analyse). Furthermore, we have to define what energy | |
| 584 | -range is covered, and what time interval is spanned by the data. All | |
| 585 | -this is done by the executable <span style="font-family: Courier New,Courier,monospace;">ctselect</span>, | |
| 586 | -which replaces the <span style="font-family: Courier New,Courier,monospace;">ctbin</span> | |
| 587 | -step in an unbinned analysis.</p> | |
| 588 | -<br> | |
| 589 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><span style="font-family: Courier New,Courier,monospace;">ctselect</span> | |
| 590 | -performs an event selection by choosing only events within a | |
| 591 | -given | |
| 592 | -region-of-interest (ROI), within a given energy band, and within a | |
| 593 | -given time interval from the input event list. The ROI is a circular | |
| 594 | -region on the sky, for which we define the centre (in celestial | |
| 595 | -coordinates) and the radius. Such a circular ROI is sometimes also | |
| 596 | -called an <span style="font-weight: bold;">acceptance cone</span>. | |
| 597 | -The following example shows how <span style="font-family: Courier New,Courier,monospace;">ctselect</span> | |
| 598 | -is invoked on the simulated | |
| 599 | -test data: </p> | |
| 600 | -<p class="Code"> | |
| 601 | -> ctselect<br> | |
| 602 | -Input event data file [test/data/crab_events.fits.gz] | |
| 603 | -events.fits<br> | |
| 604 | -Output event data file [selected_events.fits] <br> | |
| 605 | -RA for ROI centre (degrees) (0-360) [83.63] <br> | |
| 606 | -Dec for ROI centre (degrees) (-90-90) [22.01] <br> | |
| 607 | -Radius of ROI (degrees) (0-180) [3.0] <br> | |
| 608 | -Start time (MET in s) (0) [0.0] <br> | |
| 609 | -End time (MET in s) (0) [1800.0] <br> | |
| 610 | -Lower energy limit (TeV) (0) [0.1] <br> | |
| 611 | -Upper energy limit (TeV) (0) [100.0]<br> | |
| 612 | -</p> | |
| 613 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><span style="font-family: Courier New,Courier,monospace;">ctselect</span> | |
| 614 | -takes the input event list <span style="font-family: Courier New,Courier,monospace;">events.fits</span>, | |
| 615 | -performs an event selection, and writes the selected event into the | |
| 616 | -file <span style="font-family: Courier New,Courier,monospace;">selected_events.fits</span>. | |
| 617 | -The parameters it will query for are the centre of the ROI, the | |
| 618 | -radius of the ROI, the start and stop time (in seconds), and the energy | |
| 619 | -lower and upper limits (in TeV). The event selection information is | |
| 620 | -also written as a set of <span style="font-weight: bold;">data | |
| 621 | -selection keywords</span> to the output events file <span style="font-family: Courier New,Courier,monospace;">selected_events.fits</span>, | |
| 622 | -by respecting the same syntax that has been implemented for Fermi-LAT. | |
| 623 | -The following image is a screen dump of the data | |
| 624 | -selection keywords that have been written to the <span style="font-family: Courier New,Courier,monospace;">EVENTS</span> | |
| 625 | -header in the file <span style="font-family: Courier New,Courier,monospace;">selected_events.fits</span>:</p> | |
| 626 | -<br> | |
| 627 | -<p style="padding-top: 0pt; text-align: center;" class="Body"> <img style="width: 456px; height: 310px;" alt="Data selection keywords" src="images/dskeys.jpg"> | |
| 628 | -</p> | |
| 629 | -<br> | |
| 630 | -<p style="padding-top: 0pt; text-align: left;" class="Body">It | |
| 631 | -is important for unbinned analysis that these data selection keywords | |
| 632 | -exist. If they don't exist, <span style="font-family: Courier New,Courier,monospace;">ctlike</span> | |
| 633 | -will not execute in unbinned mode.</p> | |
| 634 | -<br> | |
| 635 | -<p style="padding-top: 0pt; text-align: left;" class="Body">Below | |
| 636 | -we show also for completeness the last few lines of the <span style="font-family: Courier New,Courier,monospace;">ctselect.log</span> | |
| 637 | -file. It gives information about the specific data selection that has | |
| 638 | -been performed. | |
| 639 | -</p> | |
| 640 | -<p class="Code"> | |
| 641 | -2011-02-01T19:58:34: +=================+<br> | |
| 642 | -2011-02-01T19:58:34: | Event selection |<br> | |
| 643 | -2011-02-01T19:58:34: +=================+<br> | |
| 644 | -2011-02-01T19:58:34: Time range ................: 0-1800<br> | |
| 645 | -2011-02-01T19:58:34: Energy range ..............: 0.1-100 TeV<br> | |
| 646 | -2011-02-01T19:58:34: Acceptance cone centre ....: RA=83.63, | |
| 647 | -DEC=22.01 deg<br> | |
| 648 | -2011-02-01T19:58:34: Acceptance cone radius ....: 3 deg<br> | |
| 649 | -2011-02-01T19:58:34: cfitsio selection .........: TIME | |
| 650 | ->= 0 && TIME <= 1800 && | |
| 651 | -ENERGY >= 0.1 && ENERGY <= 100 | |
| 652 | -&& ANGSEP(83.63,<br> | |
| 653 | -22.01,RA,DEC) <= 3<br> | |
| 654 | -2011-02-01T19:58:34: FITS filename .............: | |
| 655 | -events.fits[EVENTS][TIME >= 0 && TIME <= | |
| 656 | -1800 && ENERGY >= 0.1 && ENERGY | |
| 657 | -<= <br> | |
| 658 | -100 && ANGSEP(83.63,22.01,RA,DEC) <= 3]<br> | |
| 659 | -2011-02-01T19:58:34: <br> | |
| 660 | -2011-02-01T19:58:34: Application "ctselect" terminated after consuming | |
| 661 | -10 seconds of CPU time.<br> | |
| 662 | - <br> | |
| 663 | -</p> | |
| 664 | -<p style="padding-top: 0pt; text-align: left;" class="Body">Now | |
| 665 | -that the event selection has been performed, <span style="font-family: Courier New,Courier,monospace;">ctlike</span> | |
| 666 | -can be run in unbinned mode. To do this, you have to type <span style="font-family: Courier New,Courier,monospace;">UNBINNED</span> | |
| 667 | -when the maximum likelihood method is queried. The remaining parameters | |
| 668 | -are similar to binned mode, except for the last one that asks now for | |
| 669 | -an event data file instead of a counts map. Specify here the selected | |
| 670 | -event list <span style="font-family: Courier New,Courier,monospace;">selected_events.fits</span> | |
| 671 | -that has been produced by <span style="font-family: Courier New,Courier,monospace;">ctselect</span>. | |
| 672 | -To run <span style="font-family: Courier New,Courier,monospace;">ctlike</span> | |
| 673 | -in unbinned mode, type the following: | |
| 674 | -</p> | |
| 675 | -<p class="Code"> | |
| 676 | -> ctlike<br> | |
| 677 | -Event list, counts map or observation definition | |
| 678 | -file [test/data/crab_events.fits.gz] selected_events.fits<br> | |
| 679 | -Calibration database [test/irf] /usr/local/gamma/share/caldb/cta<br> | |
| 680 | -Instrument response function [kb_E_50h_v3] <br> | |
| 681 | -Source model [test/data/crab.xml] /usr/local/gamma/share/models/crab.xml<br> | |
| 682 | -Source model output file [crab_results.xml] </p> | |
| 683 | -<p style="padding-top: 0pt; text-align: left;" class="Body">You | |
| 684 | -will recognise that <span style="font-family: Courier New,Courier,monospace;">ctlike</span> | |
| 685 | -runs much faster in unbinned mode compared to binned mode. This is | |
| 686 | -understandable as the selected event list contains only 4048 events, | |
| 687 | -while the binned counts map we used before had 200 x 200 x 20 = 800000 | |
| 688 | -pixels. As unbinned maximum likelihood fitting loops over the events | |
| 689 | -(while binned maximum likelihood loops over the pixels), there are much | |
| 690 | -less operations to perform in unbinned than in binned mode (there is | |
| 691 | -some additional overhead in unbinned mode that comes from integrating | |
| 692 | -the models over the region of interest, yet this is negligible compared | |
| 693 | -to the operations needed when looping over the pixels). So as long as | |
| 694 | -you work with | |
| 695 | -short event lists, unbinned mode is faster. Unbinned <span style="font-family: Courier New,Courier,monospace;">ctlike</span> | |
| 696 | -should also be more | |
| 697 | -precise as no binning is performed, hence there is no loss of | |
| 698 | -information due to histogramming.</p> | |
| 699 | -<br> | |
| 700 | -<p style="padding-top: 0pt; text-align: left;" class="Body">Below | |
| 701 | -you see the corresponding output from the <span style="font-family: Courier New,Courier,monospace;">ctlike.log</span> | |
| 702 | -file. The fitted parameters are essentially identical to the ones found | |
| 703 | -in binned mode. The slight | |
| 704 | -difference with respect to the binned analysis may be explained by the | |
| 705 | -different event sample that has been used for the analysis: while | |
| 706 | -binned likelihood works on rectangular counts maps, unbinned likelihood | |
| 707 | -works on circular event selection regions. It is thus not possible to | |
| 708 | -select exactly the same events for both analyses!</p> | |
| 709 | -<p class="Code"> | |
| 710 | -2011-02-02T14:32:54: +=================================+<br> | |
| 711 | -2011-02-02T14:32:54: | Maximum likelihood optimisation |<br> | |
| 712 | -2011-02-02T14:32:54: +=================================+<br> | |
| 713 | -2011-02-02T14:32:54: Initial iteration: func=33730.9, | |
| 714 | -Lambda=0.001<br> | |
| 715 | -2011-02-02T14:32:54: Iteration 1: func=33728.3, | |
| 716 | -Lambda=0.0001, delta=2.56377<br> | |
| 717 | -2011-02-02T14:32:54: Iteration 2: func=33728.3, Lambda=1e-05, | |
| 718 | -delta=0.00443209<br> | |
| 719 | -2011-02-02T14:32:54: Iteration 3: func=33728.3, Lambda=1e-06, | |
| 720 | -delta=5.81336e-06<br> | |
| 721 | -2011-02-02T14:32:54: Iteration 4: func=33728.3, Lambda=1e-07, | |
| 722 | -delta=1.31549e-08<br> | |
| 723 | -2011-02-02T14:32:54: <br> | |
| 724 | -2011-02-02T14:32:54: +=========================================+<br> | |
| 725 | -2011-02-02T14:32:54: | Maximum likelihood optimisation results |<br> | |
| 726 | -2011-02-02T14:32:54: +=========================================+<br> | |
| 727 | -2011-02-02T14:32:54: === GOptimizerLM ===<br> | |
| 728 | -2011-02-02T14:32:54: Optimized function value ..: 33728.3<br> | |
| 729 | -2011-02-02T14:32:54: Absolute precision ........: 1e-06<br> | |
| 730 | -2011-02-02T14:32:54: Optimization status .......: converged<br> | |
| 731 | -2011-02-02T14:32:54: Number of parameters ......: 11<br> | |
| 732 | -2011-02-02T14:32:54: Number of free parameters .: 5<br> | |
| 733 | -2011-02-02T14:32:54: Number of iterations ......: 4<br> | |
| 734 | -2011-02-02T14:32:54: Lambda ....................: 1e-07<br> | |
| 735 | -2011-02-02T14:32:54: Maximum log likelihood ....: -33728.3<br> | |
| 736 | -2011-02-02T14:32:54: Npred .....................: 4048<br> | |
| 737 | -2011-02-02T14:32:54: === GModels ===<br> | |
| 738 | -2011-02-02T14:32:54: Number of models ..........: 2<br> | |
| 739 | -2011-02-02T14:32:54: Number of parameters ......: 11<br> | |
| 740 | -2011-02-02T14:32:54: === GModelPointSource ===<br> | |
| 741 | -2011-02-02T14:32:54: Name ......................: MAGIC Crab | |
| 742 | -(Albert et al. 2008, ApJ, 674, 1037)<br> | |
| 743 | -2011-02-02T14:32:54: Instruments ...............: all<br> | |
| 744 | -2011-02-02T14:32:54: Number of parameters ......: 6<br> | |
| 745 | -2011-02-02T14:32:54: Number of spatial par's ...: 2<br> | |
| 746 | -2011-02-02T14:32:54: RA .......................: | |
| 747 | -83.6331 [-360,360] deg (fixed,scale=1)<br> | |
| 748 | -2011-02-02T14:32:54: DEC ......................: | |
| 749 | -22.0145 [-90,90] deg (fixed,scale=1)<br> | |
| 750 | -2011-02-02T14:32:54: Number of spectral par's ..: 3<br> | |
| 751 | -2011-02-02T14:32:54: Prefactor ................: | |
| 752 | -5.83432e-16 +/- 2.44527e-17 [1e-23,1e-13] ph/cm2/s/MeV | |
| 753 | -(free,scale=1e-16)<br> | |
| 754 | -2011-02-02T14:32:54: Index ....................: | |
| 755 | --2.50725 +/- 0.0328297 [-0,-5] (free,scale=-1)<br> | |
| 756 | -2011-02-02T14:32:54: PivotEnergy ..............: | |
| 757 | -300000 [10000,1e+09] MeV (fixed,scale=1e+06)<br> | |
| 758 | -2011-02-02T14:32:54: Number of temporal par's ..: 1<br> | |
| 759 | -2011-02-02T14:32:54: Constant .................: 1 | |
| 760 | -(relative value) (fixed,scale=1)<br> | |
| 761 | -2011-02-02T14:32:54: === GCTAModelRadialAcceptance ===<br> | |
| 762 | -2011-02-02T14:32:54: Name ......................: Background<br> | |
| 763 | -2011-02-02T14:32:54: Instruments ...............: CTA<br> | |
| 764 | -2011-02-02T14:32:54: Model type ................: "Gaussian" | |
| 765 | -* "PowerLaw" * "Constant"<br> | |
| 766 | -2011-02-02T14:32:54: Number of parameters ......: 5<br> | |
| 767 | -2011-02-02T14:32:54: Number of radial par's ....: 1<br> | |
| 768 | -2011-02-02T14:32:54: Sigma ....................: | |
| 769 | -3.03472 +/- 0.0398509 [0.01,10] deg2 (free,scale=1)<br> | |
| 770 | -2011-02-02T14:32:54: Number of spectral par's ..: 3<br> | |
| 771 | -2011-02-02T14:32:54: Prefactor ................: | |
| 772 | -6.18897e-05 +/- 1.74306e-06 [0,0.001] ph/cm2/s/MeV (free,scale=1e-06)<br> | |
| 773 | -2011-02-02T14:32:54: Index ....................: | |
| 774 | --1.86369 +/- 0.0156024 [-0,-5] (free,scale=-1)<br> | |
| 775 | -2011-02-02T14:32:54: PivotEnergy ..............: | |
| 776 | -1e+06 [10000,1e+09] MeV (fixed,scale=1e+06)<br> | |
| 777 | -2011-02-02T14:32:54: Number of temporal par's ..: 1<br> | |
| 778 | -2011-02-02T14:32:54: Constant .................: 1 | |
| 779 | -(relative value) (fixed,scale=1)<br> | |
| 780 | -2011-02-02T14:32:54: <br> | |
| 781 | -2011-02-02T14:32:54: Application "ctlike" terminated after consuming 2 | |
| 782 | -seconds of CPU time.<br> | |
| 783 | -</p> | |
| 784 | -<br> | |
| 785 | -<br> | |
| 786 | -<p class="Title"><a name="Some_final_remarks"></a>Some | |
| 787 | -final remarks <span style="color: red;">(updated for | |
| 788 | -ctools-00-02-03)</span></p> | |
| 789 | -<p style="padding-top: 0pt; text-align: left;" class="Body">First | |
| 790 | -you should notice that <span style="font-weight: bold;">ctools | |
| 791 | -are in the prototyping phase</span>. If GammaLib and ctatools are | |
| 792 | -installed on your system you should be able to run the above example | |
| 793 | -without any problems, and you can do even much more.</p> | |
| 794 | -<br> | |
| 795 | -<p style="padding-top: 0pt; text-align: left;" class="Body">You | |
| 796 | -may specify longer integration times for the simulation and for the | |
| 797 | -analysis, you can have several point sources in the | |
| 798 | -XML file, and you may try to fit all sources simultaneously using <span style="font-family: Courier New,Courier,monospace;">ctlike</span>. | |
| 799 | -Yet, for the moment, | |
| 800 | -you are limited to point sources (<span style="font-family: Courier New,Courier,monospace;">type="PointSource"</span>) | |
| 801 | -and their spectra should be power laws (<span style="font-family: Courier New,Courier,monospace;">type="PowerLaw"</span>). | |
| 802 | -Although there are already a number of other spectral types available | |
| 803 | -in GammaLib, I have not yet had the time to implement their Monte Carlo | |
| 804 | -generator, hence you cannot simulate events for spectral laws different | |
| 805 | -from power laws for the moment. You may try, however, fitting the | |
| 806 | -sources with other spectral laws, yet I have not tried this, and I | |
| 807 | -cannot guarantee that it works.</p> | |
| 808 | -<br> | |
| 809 | -<p style="padding-top: 0pt; text-align: left;" class="Body">There | |
| 810 | -is also only one radial acceptance model defined so far, and you should | |
| 811 | -use the model definition that is given in the example.</p> | |
| 812 | -<br> | |
| 813 | -<p style="padding-top: 0pt; text-align: left;" class="Body">So | |
| 814 | -far, there are four projection types supported: <span style="font-family: Courier New,Courier,monospace;">AZP</span>, | |
| 815 | -<span style="font-family: Courier New,Courier,monospace;">CAR</span>, | |
| 816 | -<span style="font-family: Courier New,Courier,monospace;">STG</span>, | |
| 817 | -and <span style="font-family: Courier New,Courier,monospace;">TAN</span>. | |
| 818 | -You can use them in Galactic (<span style="font-family: Courier New,Courier,monospace;">GAL</span>) | |
| 819 | -or celestial (<span style="font-family: Courier New,Courier,monospace;">CEL</span>) | |
| 820 | -coordinates.</p> | |
| 821 | -<br> | |
| 822 | -<p style="padding-top: 0pt; text-align: left;" class="Body">The | |
| 823 | -instrument response function that has been implemented so far in | |
| 824 | -GammaLib (yes, it is in GammaLib, not in ctatools) is very simple. It | |
| 825 | -just reads Konrad's effective area files to get the onaxis effective | |
| 826 | -area. It is assumed that the effective area varies with offset angle, | |
| 827 | -and a Gaussian as function of offset angle squared has been implemented | |
| 828 | -for the functional form of this variation. By default, the <span style="font-family: Courier New,Courier,monospace;">sigma</span> | |
| 829 | -of this Gaussian is fixed to 3 deg<sup>2</sup>. | |
| 830 | -This value can not be modified using the ctatools, but it can be | |
| 831 | -altered using the Python interface. In any case, this treatment of the | |
| 832 | -offset angle dependence is only preliminary and <span style="font-weight: bold;">will be changed as soon as a | |
| 833 | -format has been defined for the CTA response functions</span>.</p> | |
| 834 | -<p style="padding-top: 0pt; text-align: left;" class="Body"></p> | |
| 835 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><span style="font-family: Courier New,Courier,monospace;">ctobssim</span> | |
| 836 | -and <span style="font-family: Courier New,Courier,monospace;">ctlike</span> | |
| 837 | -query for the name of the instrument response function file, and you | |
| 838 | -may specify different | |
| 839 | -response files if you want to study the impact of array configuration | |
| 840 | -on the results. To see what is available, type (adjust the | |
| 841 | -path corresponding to your installation directory): | |
| 842 | -</p> | |
| 843 | -<p class="Code">> ls | |
| 844 | -/usr/local/gamma/share/caldb/cta<br> | |
| 845 | -kb_A_50h_v3.dat | |
| 846 | -kb_D_50h_v3.dat | |
| 847 | -kb_G_50h_v3.dat | |
| 848 | -kb_J_50h_v3.dat | |
| 849 | -kb_NB_50h_v3.dat<br> | |
| 850 | -kb_B_50h_v3.dat | |
| 851 | -kb_E_50h_v3.dat | |
| 852 | -kb_H_50h_v3.dat | |
| 853 | -kb_K_50h_v3.dat<br> | |
| 854 | -kb_C_50h_v3.dat | |
| 855 | -kb_F_50h_v3.dat | |
| 856 | -kb_I_50h_v3.dat | |
| 857 | -kb_NA_50h_v3.dat<br> | |
| 858 | -</p> | |
| 859 | -<p style="padding-top: 0pt; text-align: left;" class="Body">The | |
| 860 | -point spread function (PSF) definition is also read from these files. | |
| 861 | -So far, the <span style="font-weight: bold;">PSF is a | |
| 862 | -simple Gaussian</span>. The energy dependent Gaussian | |
| 863 | -width is defined by | |
| 864 | -the r68 parameter in Konrad's response files.</p> | |
| 865 | -<br> | |
| 866 | -<p style="padding-top: 0pt; text-align: left;" class="Body"><span style="font-weight: bold;">Lifetime | |
| 867 | -is assumed to be 100%</span>. I'm not sure whether lifetime is | |
| 868 | -folded in the | |
| 869 | -effective areas of Konrad's response files, if not, a corresponding | |
| 870 | -correction should be applied to all source fluxes obtained with the | |
| 871 | -ctatools. | |
| 872 | -</p> | |
| 873 | -<br> | |
| 874 | -<p style="padding-top: 0pt; text-align: left;" class="Body">Finally, | |
| 875 | -if | |
| 876 | -you want to have the log file also dumped to the screen as the | |
| 877 | -executable runs, add <span style="font-family: Courier New,Courier,monospace;">debug=yes</span> | |
| 878 | -on the command line, e.g. | |
| 879 | -</p> | |
| 880 | -<p class="Code">> ctlike debug=yes | |
| 881 | -</p> | |
| 882 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 883 | -And now: <span style="font-weight: bold;">enjoy playing | |
| 884 | -with ctatools</span>. | |
| 885 | -</p> | |
| 886 | -<br> | |
| 887 | -<br> | |
| 888 | -<p class="Title"><a name="Getting_Help"></a>Getting | |
| 889 | -Help</p> | |
| 890 | -<p style="padding-top: 0pt; text-align: left;" class="Body"> | |
| 891 | -Any questions, bug reports, or suggested enhancements related to | |
| 892 | -ctatools and GammaLib should be submitted | |
| 893 | -to <a href="mailto:jurgen.knodlseder@irap.omp.eu">jurgen.knodlseder@irap.omp.eu</a>. | |
| 894 | -</p> | |
| 895 | -<br> | |
| 896 | -<hr size="1"> | |
| 897 | -<p class="Body">Last updated by: Jürgen | |
| 898 | -Knödlseder (4 April 2012) | |
| 899 | -</p> | |
| 900 | -</body></html> | |
| 901 | 0 | \ No newline at end of file |
doc/html/images/example_survey.jpg renamed to doc/source/_static/example_survey.jpg
245 KB
test/example_survey.png renamed to doc/source/_static/example_survey.png
251 KB
doc/html/images/header.jpg renamed to doc/source/_static/header.jpg
35.5 KB
doc/html/images/survey-pop2.jpg renamed to doc/source/_static/survey-pop2.jpg
350 KB
doc/source/_templates/index.html
| ... | ... | @@ -37,7 +37,7 @@ CTA format), Fermi/LAT and COMPTEL data. {%endtrans%}</p> |
| 37 | 37 | <li><a class="smalllink" href="{{ pathto("user_manual/getting_started/quickstart") }}">{%trans%}Quickstart{%endtrans%}</a></li> |
| 38 | 38 | <li><a class="smalllink" href="{{ pathto("user_manual/getting_started/response") }}">{%trans%}CTA response functions{%endtrans%}</a></li> |
| 39 | 39 | <li><a class="smalllink" href="{{ pathto("user_manual/index") }}">{%trans%}User Manual{%endtrans%}</a></li> |
| 40 | -<li><a class="smalllink" href="{{ pathto("documentation") }}">{%trans%}ctools reference{%endtrans%}</a></li> | |
| 40 | +<li><a class="smalllink" href="{{ pathto("tools_reference") }}">{%trans%}ctools Reference{%endtrans%}</a></li> | |
| 41 | 41 | <li><a class="smalllink" href="doxygen/index.html">{%trans%}Source code documentation (Doxygen){%endtrans%}</a></li> |
| 42 | 42 | </ul> |
| 43 | 43 | ... | ... |
doc/source/conf.py
| ... | ... | @@ -23,9 +23,34 @@ import sys, os |
| 23 | 23 | # If your documentation needs a minimal Sphinx version, state it here. |
| 24 | 24 | #needs_sphinx = '1.0' |
| 25 | 25 | |
| 26 | +extensions = ['sphinx.ext.intersphinx', 'sphinx.ext.mathjax'] | |
| 27 | + | |
| 26 | 28 | # Add any Sphinx extension module names here, as strings. They can be extensions |
| 27 | 29 | # coming with Sphinx (named 'sphinx.ext.*') or your custom ones. |
| 28 | -extensions = ['sphinx.ext.mathjax'] | |
| 30 | +try: | |
| 31 | + import sphinxcontrib.doxylink | |
| 32 | + extensions += ['sphinxcontrib.doxylink'] | |
| 33 | +except: | |
| 34 | + print('WARNING: Python package `sphinxcontrib.doxylink` not found.\n' | |
| 35 | + 'WARNING: Please install it to get fully-functional Sphinx docs.\n' | |
| 36 | + 'WARNING: https://pypi.python.org/pypi/sphinxcontrib-doxylink\n') | |
| 37 | + | |
| 38 | +# http://sphinx-doc.org/extensions.html#where-to-put-your-own-extensions | |
| 39 | +sys.path.append(os.path.abspath('exts')) | |
| 40 | + | |
| 41 | +# We bundled the `rawfiles` sphinx extension to copy over files on sphinx build | |
| 42 | +# https://pypi.python.org/pypi/sphinxcontrib-rawfiles | |
| 43 | +extensions += ['rawfiles'] | |
| 44 | + | |
| 45 | +# You can list any files or folders here that should be copied 1:1 to the html | |
| 46 | +# output directory by sphinx-build | |
| 47 | +rawfiles = ['fhelp'] | |
| 48 | + | |
| 49 | +# This makes it possible to cross-link to the GammaLib Sphinx docs | |
| 50 | +# see http://sphinx-doc.org/ext/intersphinx.html | |
| 51 | +intersphinx_mapping = dict() | |
| 52 | +intersphinx_mapping['gammalib'] = ('http://gammalib.sourceforge.net/', None) | |
| 53 | + | |
| 29 | 54 | |
| 30 | 55 | # Add any paths that contain templates here, relative to this directory. |
| 31 | 56 | templates_path = ['_templates'] |
| ... | ... | @@ -241,3 +266,9 @@ texinfo_documents = [ |
| 241 | 266 | |
| 242 | 267 | # How to display URL addresses: 'footnote', 'no', or 'inline'. |
| 243 | 268 | #texinfo_show_urls = 'footnote' |
| 269 | + | |
| 270 | +# Define connection to Doxygen API docs | |
| 271 | +# http://pythonhosted.org/sphinxcontrib-doxylink/index.html#confval-doxylink | |
| 272 | +doxylink = { | |
| 273 | + 'doxy' : ('doc/html/doxygen/ctools.tag', 'doxygen/'), | |
| 274 | +} | ... | ... |
doc/source/exts/rawfiles.py
0 → 100644
| 1 | +# -*- coding: utf-8 -*- | |
| 2 | +# Copied from https://bitbucket.org/birkenfeld/sphinx-contrib/ | |
| 3 | +import os | |
| 4 | +import shutil | |
| 5 | + | |
| 6 | +def on_html_collect_pages(app): | |
| 7 | + for f in app.builder.config.rawfiles: | |
| 8 | + src = os.path.join(app.srcdir, f) | |
| 9 | + dst = os.path.join(app.builder.outdir, f) | |
| 10 | + # cleanup | |
| 11 | + if os.path.exists(dst): | |
| 12 | + if os.path.isdir(dst): | |
| 13 | + shutil.rmtree(dst) | |
| 14 | + elif os.path.isfile(dst): | |
| 15 | + os.remove(dst) | |
| 16 | + # copy file | |
| 17 | + if os.path.exists(src): | |
| 18 | + if os.path.isdir(src): | |
| 19 | + shutil.copytree(src, dst) | |
| 20 | + elif os.path.isfile(src): | |
| 21 | + shutil.copy(src, dst) | |
| 22 | + else: | |
| 23 | + msg = ('rawfiles target is not directory or file "%s".' % src) | |
| 24 | + app.builder.warn(msg) | |
| 25 | + else: | |
| 26 | + msg = ('rawfiles can not find "%s".' % src) | |
| 27 | + app.builder.warn(msg) | |
| 28 | + | |
| 29 | + return () | |
| 30 | + | |
| 31 | +def setup(app): | |
| 32 | + app.add_config_value('rawfiles', [], 'html') | |
| 33 | + app.connect("html-collect-pages", on_html_collect_pages) | |
| 0 | 34 | \ No newline at end of file | ... | ... |
doc/html/doc/cspull.html renamed to doc/source/fhelp/cspull.html
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | <html lang="en"><head> |
| 3 | 3 | <meta content="text/html; charset=ISO-8859-1" http-equiv="content-type"><title>cspull</title> |
| 4 | 4 | |
| 5 | -<link rel="stylesheet" type="text/css" media="screen,print" href="../ctools.css"> | |
| 5 | +<link rel="stylesheet" type="text/css" media="screen,print" href="ctools.css"> | |
| 6 | 6 | <meta content="Jurgen Knodlseder" name="author"> |
| 7 | 7 | <meta content="CTA pull distribution generator" name="description"> |
| 8 | 8 | </head> |
| ... | ... | @@ -16,7 +16,7 @@ parameters.</p> |
| 16 | 16 | pull distributions for all parameters in a model. The pull is defined |
| 17 | 17 | by the fitted value minus the simulated true value, devided by the |
| 18 | 18 | statistical uncertainty. For unbiased parameter estimates, the pull |
| 19 | -distribution should follow a Gaussian centred on 0 and a sigma | |
| 19 | +distribution should follow a Gaussian centred on 0 and a sigma | |
| 20 | 20 | parameter of 1.</p> |
| 21 | 21 | <br> |
| 22 | 22 | <p class="Title">Links</p> | ... | ... |
doc/html/doc/cssens.html renamed to doc/source/fhelp/cssens.html
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | <html lang="en"><head> |
| 3 | 3 | <meta content="text/html; charset=ISO-8859-1" http-equiv="content-type"><title>cssens</title> |
| 4 | 4 | |
| 5 | -<link rel="stylesheet" type="text/css" media="screen,print" href="../ctools.css"> | |
| 5 | +<link rel="stylesheet" type="text/css" media="screen,print" href="ctools.css"> | |
| 6 | 6 | <meta content="Jurgen Knodlseder" name="author"> |
| 7 | 7 | <meta content="CTA sensitivity estimator" name="description"> |
| 8 | 8 | </head> | ... | ... |
doc/html/doc/cstsdist.html renamed to doc/source/fhelp/cstsdist.html
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | <html lang="en"><head> |
| 3 | 3 | <meta content="text/html; charset=ISO-8859-1" http-equiv="content-type"><title>cstsdist</title> |
| 4 | 4 | |
| 5 | -<link rel="stylesheet" type="text/css" media="screen,print" href="../ctools.css"> | |
| 5 | +<link rel="stylesheet" type="text/css" media="screen,print" href="ctools.css"> | |
| 6 | 6 | <meta content="Jurgen Knodlseder" name="author"> |
| 7 | 7 | <meta content="CTA TS distribution generator" name="description"> |
| 8 | 8 | </head> | ... | ... |
doc/html/doc/ctbin.html renamed to doc/source/fhelp/ctbin.html
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | <html lang="en"><head> |
| 3 | 3 | <meta content="text/html; charset=ISO-8859-1" http-equiv="content-type"><title>ctbin</title> |
| 4 | 4 | |
| 5 | -<link rel="stylesheet" type="text/css" media="screen,print" href="../ctools.css"> | |
| 5 | +<link rel="stylesheet" type="text/css" media="screen,print" href="ctools.css"> | |
| 6 | 6 | <meta content="Jurgen Knodlseder" name="author"> |
| 7 | 7 | <meta content="CTA event binning tool" name="description"> |
| 8 | 8 | </head> | ... | ... |
doc/html/doc/ctlike.html renamed to doc/source/fhelp/ctlike.html
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | <html lang="en"><head> |
| 3 | 3 | <meta content="text/html; charset=ISO-8859-1" http-equiv="content-type"><title>ctlike</title> |
| 4 | 4 | |
| 5 | -<link rel="stylesheet" type="text/css" media="screen,print" href="../ctools.css"> | |
| 5 | +<link rel="stylesheet" type="text/css" media="screen,print" href="ctools.css"> | |
| 6 | 6 | <meta content="Jurgen Knodlseder" name="author"> |
| 7 | 7 | <meta content="CTA maximum likelihood fitting tool" name="description"></head> |
| 8 | 8 | <body> | ... | ... |
doc/html/doc/ctmodel.html renamed to doc/source/fhelp/ctmodel.html
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | <html lang="en"><head> |
| 3 | 3 | <meta content="text/html; charset=ISO-8859-1" http-equiv="content-type"><title>ctmodel</title> |
| 4 | 4 | |
| 5 | -<link rel="stylesheet" type="text/css" media="screen,print" href="../ctools.css"> | |
| 5 | +<link rel="stylesheet" type="text/css" media="screen,print" href="ctools.css"> | |
| 6 | 6 | <meta content="Jurgen Knodlseder" name="author"> |
| 7 | 7 | <meta content="CTA modelling tool" name="description"> |
| 8 | 8 | </head> | ... | ... |
doc/html/doc/ctobssim.html renamed to doc/source/fhelp/ctobssim.html
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | <html lang="en"><head> |
| 3 | 3 | <meta content="text/html; charset=ISO-8859-1" http-equiv="content-type"><title>ctobssim</title> |
| 4 | 4 | |
| 5 | -<link rel="stylesheet" type="text/css" media="screen,print" href="../ctools.css"> | |
| 5 | +<link rel="stylesheet" type="text/css" media="screen,print" href="ctools.css"> | |
| 6 | 6 | <meta content="Jurgen Knodlseder" name="author"> |
| 7 | 7 | <meta content="CTA observation simulation tool" name="description"> |
| 8 | 8 | </head> | ... | ... |
doc/html/ctools.css renamed to doc/source/fhelp/ctools.css
doc/html/doc/ctselect.html renamed to doc/source/fhelp/ctselect.html
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | <html lang="en"><head> |
| 3 | 3 | <meta content="text/html; charset=ISO-8859-1" http-equiv="content-type"><title>ctselect</title> |
| 4 | 4 | |
| 5 | -<link rel="stylesheet" type="text/css" media="screen,print" href="../ctools.css"> | |
| 5 | +<link rel="stylesheet" type="text/css" media="screen,print" href="ctools.css"> | |
| 6 | 6 | <meta content="Jurgen Knodlseder" name="author"> |
| 7 | 7 | <meta content="CTA event sub-selection tool" name="description"> |
| 8 | 8 | </head> | ... | ... |
doc/html/doc/ctskymap.html renamed to doc/source/fhelp/ctskymap.html
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | <html lang="en"><head> |
| 3 | 3 | <meta content="text/html; charset=ISO-8859-1" http-equiv="content-type"><title>ctskymap</title> |
| 4 | 4 | |
| 5 | -<link rel="stylesheet" type="text/css" media="screen,print" href="../ctools.css"> | |
| 5 | +<link rel="stylesheet" type="text/css" media="screen,print" href="ctools.css"> | |
| 6 | 6 | <meta content="Jurgen Knodlseder" name="author"> |
| 7 | 7 | <meta content="CTA sky mapping tool" name="description"> |
| 8 | 8 | </head> | ... | ... |
doc/html/doc/ftools_usage.html renamed to doc/source/fhelp/ftools_usage.html
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | <html lang="en"><head> |
| 3 | 3 | <meta content="text/html; charset=ISO-8859-1" http-equiv="content-type"><title>Basic FTOOL Parameter Usage</title> |
| 4 | 4 | |
| 5 | -<link rel="stylesheet" type="text/css" media="screen,print" href="../ctools.css"> | |
| 5 | +<link rel="stylesheet" type="text/css" media="screen,print" href="ctools.css"> | |
| 6 | 6 | <meta content="Jurgen Knodlseder" name="author"> |
| 7 | 7 | <meta content="Basic FTOOL Parameter Usage" name="description"> |
| 8 | 8 | </head> | ... | ... |
doc/source/tools_reference.rst
0 → 100644
| 1 | +ctools reference manual | |
| 2 | +======================= | |
| 3 | + | |
| 4 | +General information on FTOOLS usage can be found `here <fhelp/ftools_usage.html>`__. | |
| 5 | + | |
| 6 | +Below you can find links to the command line tools help pages that are available in ctools. | |
| 7 | + | |
| 8 | +ctools | |
| 9 | +------ | |
| 10 | + | |
| 11 | +* `ctbin <fhelp/ctbin.html>`__ --- Calculates counts maps for binned maximum likelihood analysis. | |
| 12 | +* `ctlike <fhelp/ctlike.html>`__ --- Performs binned and unbinned maximum likelihood analysis of CTA data. | |
| 13 | +* `ctmodel <fhelp/ctmodel.html>`__ --- Calculates model counts maps. | |
| 14 | +* `ctobssim <fhelp/ctobssim.html>`__ --- Observation simulator. | |
| 15 | +* `ctselect <fhelp/ctselect.html>`__ --- Performs sub-selections of event data which are needed for unbinned maximum likelihood analysis. | |
| 16 | +* `ctskymap <fhelp/ctskymap.html>`__ --- Computes sky map from CTA data. | |
| 17 | + | |
| 18 | +cscripts | |
| 19 | +-------- | |
| 20 | + | |
| 21 | +* `cspull <fhelp/cspull.html>`__ --- Generates pull distribution for all model parameters. | |
| 22 | +* `cssens <fhelp/cssens.html>`__ --- Computes the CTA sensitivity for a number of energy bins. | |
| 23 | +* `cstsdist <fhelp/cstsdist.html>`__ --- Generates the TS distribution for a particular model. | ... | ... |
doc/source/user_manual/enhancing.rst
| 1 | 1 | Enhancing ctools |
| 2 | 2 | ================ |
| 3 | 3 | |
| 4 | +Most ctools users will run the existing ctools or cscripts, | |
| 5 | +which cover the most common simulation and analysis needs. | |
| 6 | + | |
| 7 | +Advanced users and developers can write their own tools and scripts. | |
| 8 | +The following C++ and Python examples can help get you started: | |
| 9 | + | |
| 10 | +* :doxy:`ctlike` | |
| 11 | +* :doxy:`obsutils` | |
| 12 | +* :doxy:`cstsdist` | |
| 13 | + | |
| 4 | 14 | Building your own cscripts |
| 5 | 15 | -------------------------- |
| 6 | 16 | |
| 7 | 17 | Building your own ctools |
| 8 | ------------------------- | |
| 9 | 18 | \ No newline at end of file |
| 19 | +------------------------ | ... | ... |
doc/source/user_manual/getting_started/models.rst
| ... | ... | @@ -2,4 +2,5 @@ Models available in ctools |
| 2 | 2 | -------------------------- |
| 3 | 3 | |
| 4 | 4 | To be written ... (you may read |
| 5 | -http://gammalib.sourceforge.net/user_manual/modules/model.html#spectral-models for the time being). | |
| 5 | +:ref:`the model section in the GammaLib user manual <gammalib:sec_model>` | |
| 6 | +for the time being). | ... | ... |
doc/source/user_manual/index.rst
doc/source/user_manual/reference.rst deleted
examples/make_binned_analysis.py
| ... | ... | @@ -46,11 +46,11 @@ def pipeline_v1(): |
| 46 | 46 | in a ftools-like approach. |
| 47 | 47 | """ |
| 48 | 48 | # Set script parameters |
| 49 | - model_name = "${GAMMALIB}/share/models/crab.xml" | |
| 49 | + model_name = "${CTOOLS}/share/models/crab.xml" | |
| 50 | 50 | events_name = "events.fits" |
| 51 | 51 | cntmap_name = "cntmap.fits" |
| 52 | 52 | result_name = "results.xml" |
| 53 | - caldb = "${GAMMALIB}/share/caldb/cta" | |
| 53 | + caldb = "${CALDB}/data/cta/dummy/bcf" | |
| 54 | 54 | irf = "cta_dummy_irf" |
| 55 | 55 | ra = 83.63 |
| 56 | 56 | dec = 22.01 |
| ... | ... | @@ -151,8 +151,8 @@ def pipeline_v2(): |
| 151 | 151 | At the end, results are plotted (if matplotlib is installed) |
| 152 | 152 | """ |
| 153 | 153 | # Set script parameters |
| 154 | - model_name = "${GAMMALIB}/share/models/crab.xml" | |
| 155 | - caldb = "${GAMMALIB}/share/caldb/cta" | |
| 154 | + model_name = "${CTOOLS}/share/models/crab.xml" | |
| 155 | + caldb = "${CALDB}/data/cta/dummy/bcf" | |
| 156 | 156 | irf = "cta_dummy_irf" |
| 157 | 157 | ra = 83.63 |
| 158 | 158 | dec = 22.01 | ... | ... |
examples/make_pull_at_sensitivity_limit.py
| ... | ... | @@ -61,7 +61,7 @@ def create_pull(loge, emin, emax, models, ntrials=100, duration=180000.0, \ |
| 61 | 61 | pull.models(models) |
| 62 | 62 | pull["outfile"] = outfile |
| 63 | 63 | pull["ntrials"] = ntrials |
| 64 | - pull["caldb"] = "$GAMMALIB/share/caldb/cta" | |
| 64 | + pull["caldb"] = "${CALDB}/data/cta/dummy/bcf" | |
| 65 | 65 | pull["irf"] = "cta_dummy_irf" |
| 66 | 66 | pull["ra"] = 83.6331 |
| 67 | 67 | pull["dec"] = 22.0145 | ... | ... |
examples/make_ts_distributions.py
| ... | ... | @@ -65,12 +65,12 @@ def create_ts(loge, emin, emax, ntrials=100, duration=180000.0, \ |
| 65 | 65 | tsdist = cstsdist() |
| 66 | 66 | tsdist["outfile"] = outfile |
| 67 | 67 | tsdist["ntrials"] = ntrials |
| 68 | - tsdist["caldb"] = "$GAMMALIB/share/caldb/cta" | |
| 68 | + tsdist["caldb"] = "${CALDB}/data/cta/dummy/bcf" | |
| 69 | 69 | tsdist["irf"] = "cta_dummy_irf" |
| 70 | 70 | tsdist["type"] = "point" |
| 71 | 71 | tsdist["index"] = -2.48 |
| 72 | 72 | tsdist["offset"] = 0.0 |
| 73 | - tsdist["bkg"] = "$GAMMALIB/share/models/bkg_dummy.txt" | |
| 73 | + tsdist["bkg"] = "${CTOOLS}/share/models/bkg_dummy.txt" | |
| 74 | 74 | tsdist["emin"] = float(emin) |
| 75 | 75 | tsdist["emax"] = float(emax) |
| 76 | 76 | tsdist["enumbins"] = int(enumbins) | ... | ... |
examples/make_unbinned_analysis.py
| ... | ... | @@ -49,11 +49,11 @@ def pipeline_v1(): |
| 49 | 49 | in a ftools-like approach. |
| 50 | 50 | """ |
| 51 | 51 | # Set script parameters |
| 52 | - model_name = "${GAMMALIB}/share/models/crab.xml" | |
| 52 | + model_name = "${CTOOLS}/share/models/crab.xml" | |
| 53 | 53 | events_name = "events.fits" |
| 54 | 54 | selected_events_name = "selected_events.fits" |
| 55 | 55 | result_name = "results.xml" |
| 56 | - caldb = "${GAMMALIB}/share/caldb/cta" | |
| 56 | + caldb = "${CALDB}/data/cta/dummy/bcf" | |
| 57 | 57 | irf = "cta_dummy_irf" |
| 58 | 58 | ra = 83.63 |
| 59 | 59 | dec = 22.01 |
| ... | ... | @@ -143,8 +143,8 @@ def pipeline_v2(show_data): |
| 143 | 143 | FITS files on disk. All data is only hold in memory. |
| 144 | 144 | """ |
| 145 | 145 | # Set script parameters |
| 146 | - model_name = "${GAMMALIB}/share/models/crab.xml" | |
| 147 | - caldb = "${GAMMALIB}/share/caldb/cta" | |
| 146 | + model_name = "${CTOOLS}/share/models/crab.xml" | |
| 147 | + caldb = "${CALDB}/data/cta/dummy/bcf" | |
| 148 | 148 | irf = "cta_dummy_irf" |
| 149 | 149 | ra = 83.63 |
| 150 | 150 | dec = 22.01 | ... | ... |
pyext/Makefile.am
| ... | ... | @@ -30,7 +30,7 @@ |
| 30 | 30 | if WITH_PYTHON |
| 31 | 31 | |
| 32 | 32 | # Set SWIG interface filename |
| 33 | -SWIG_INTERFACE = ctools.i | |
| 33 | +SWIG_INTERFACE = ctools/ctools.i | |
| 34 | 34 | |
| 35 | 35 | # Dependencies. Obtain the list of dependencies by executing ls. |
| 36 | 36 | DEPFILES = $(shell ls *.Pi 2> /dev/null) |
| ... | ... | @@ -39,7 +39,7 @@ include $(DEPFILES) |
| 39 | 39 | # Define the targets that should be built before anything else. |
| 40 | 40 | # These are honored by 'make all', 'make check' and 'make install'. |
| 41 | 41 | if HAVE_SWIG |
| 42 | - BUILT_SOURCES = ctools_wrap.cpp ctools.py | |
| 42 | + BUILT_SOURCES = ctools/ctools_wrap.cpp ctools/ctools.py | |
| 43 | 43 | endif |
| 44 | 44 | |
| 45 | 45 | # Set files to be removed in clean. We remove only the SWIG interface files |
| ... | ... | @@ -59,7 +59,7 @@ all-local: build |
| 59 | 59 | # of this specific directory. So we build one instance for checking, the other |
| 60 | 60 | # for installation. If we could get the name built by distutils we could avoid |
| 61 | 61 | # this. |
| 62 | -build: ctools_wrap.cpp ctools.py | |
| 62 | +build: ctools/ctools_wrap.cpp ctools/ctools.py | |
| 63 | 63 | $(PYTHON) setup.py build |
| 64 | 64 | |
| 65 | 65 | install-exec-local: |
| ... | ... | @@ -70,9 +70,9 @@ clean-local: |
| 70 | 70 | |
| 71 | 71 | # Rule for building Python wrappers using SWIG. |
| 72 | 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 | 76 | then mv -f "ctools.Tpi" "ctools.Pi"; else rm -f "ctools.Tpi"; exit 1; fi |
| 77 | 77 | endif # HAVE_SWIG |
| 78 | 78 | ... | ... |
pyext/ctexpcube.i
0 → 100644
| 1 | +/*************************************************************************** | |
| 2 | + * ctexpcube - CTA exposure cube tool * | |
| 3 | + * ----------------------------------------------------------------------- * | |
| 4 | + * copyright (C) 2014 by Juergen Knoedlseder * | |
| 5 | + * ----------------------------------------------------------------------- * | |
| 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 | +/** | |
| 22 | + * @file ctexpcube.i | |
| 23 | + * @brief CTA exposure cube tool definition | |
| 24 | + * @author Juergen Knoedlseder | |
| 25 | + */ | |
| 26 | +%{ | |
| 27 | +/* Put headers and other declarations here that are needed for compilation */ | |
| 28 | +#include "ctexpcube.hpp" | |
| 29 | +#include "GTools.hpp" | |
| 30 | +%} | |
| 31 | + | |
| 32 | + | |
| 33 | +/***********************************************************************//** | |
| 34 | + * @class ctexpcube | |
| 35 | + * | |
| 36 | + * @brief CTA exposure cube tool | |
| 37 | + ***************************************************************************/ | |
| 38 | +class ctexpcube : public GApplication { | |
| 39 | + | |
| 40 | +public: | |
| 41 | + // Constructors and destructors | |
| 42 | + ctexpcube(void); | |
| 43 | + explicit ctexpcube(GObservations obs); | |
| 44 | + ctexpcube(int argc, char *argv[]); | |
| 45 | + ctexpcube(const ctexpcube& app); | |
| 46 | + virtual ~ctexpcube(void); | |
| 47 | + | |
| 48 | + // Methods | |
| 49 | + void clear(void); | |
| 50 | + void execute(void); | |
| 51 | + void run(void); | |
| 52 | + void save(void); | |
| 53 | + GObservations& obs(void); | |
| 54 | + void get_parameters(void); | |
| 55 | + void get_obs(void); | |
| 56 | + void set_response(void); | |
| 57 | + void get_ebounds(void); | |
| 58 | + void set_from_cntmap(const std::string& filename); | |
| 59 | +}; | |
| 60 | + | |
| 61 | + | |
| 62 | +/***********************************************************************//** | |
| 63 | + * @brief CTA exposure map tool Python extension | |
| 64 | + ***************************************************************************/ | |
| 65 | +%extend ctexpcube { | |
| 66 | + ctexpcube copy() { | |
| 67 | + return (*self); | |
| 68 | + } | |
| 69 | +} | ... | ... |
pyext/ctools/__init__.py
0 → 100644
pyext/ctools.i renamed to pyext/ctools/ctools.i
| 1 | 1 | /*************************************************************************** |
| 2 | 2 | * ctools - SWIG file * |
| 3 | 3 | * ----------------------------------------------------------------------- * |
| 4 | - * copyright (C) 2010-2013 by Juergen Knoedlseder * | |
| 4 | + * copyright (C) 2010-2014 by Juergen Knoedlseder * | |
| 5 | 5 | * ----------------------------------------------------------------------- * |
| 6 | 6 | * * |
| 7 | 7 | * This program is free software: you can redistribute it and/or modify * |
| ... | ... | @@ -51,3 +51,4 @@ |
| 51 | 51 | %include "ctmodel.i" |
| 52 | 52 | %include "ctselect.i" |
| 53 | 53 | %include "ctskymap.i" |
| 54 | +%include "ctexpcube.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 | |
| 519 | + | |
| 520 | + | |
| 521 | +# ================================================= # | |
| 522 | +# Make spectrum for one Source in a CTA Observation # | |
| 523 | +# ================================================= # | |
| 524 | +def getlogbins(nbins=12,emin=0.1,emax=100): | |
| 525 | + """ | |
| 526 | + Returns list with logarithmiscenergy bins in | |
| 527 | + format [[emin1,emax1],[emin2,emax2],..] | |
| 528 | + | |
| 529 | + Parameters: | |
| 530 | + nbins Number of bins | |
| 531 | + emin Minimum energy | |
| 532 | + emax Maximum energy | |
| 533 | + """ | |
| 534 | + lgemin=numpy.log10(emin) | |
| 535 | + lgemax=numpy.log10(emax) | |
| 536 | + lgbins=[] | |
| 537 | + lgedges=numpy.arange(lgemin,lgemax,(lgemax-lgemin)/nbins) | |
| 538 | + lgedges=numpy.append(lgedges,[lgemax]) | |
| 539 | + for ii in range(len(lgedges)-1): | |
| 540 | + lgbins.append([lgedges[ii],lgedges[ii+1]]) | |
| 541 | + ebins=numpy.power(10,lgbins) | |
| 542 | + return ebins | |
| 543 | + | |
| 544 | +#Please note that specpoints can make spectra for any given | |
| 545 | +#spectral model, as long as it has a "Prefactor" parameter | |
| 546 | +#which scales the funtion | |
| 547 | +def specpoints(obs,source,ebins): | |
| 548 | + """ | |
| 549 | + Returns spectral points in the given binning for | |
| 550 | + the passed observations and source. | |
| 551 | + | |
| 552 | + Parameters: | |
| 553 | + obs - Pointing direction | |
| 554 | + source - Source name for which spectrum should be done | |
| 555 | + ebins - Binning passed as a list in units of TeV as | |
| 556 | + [[emin1,emax1],[emin2,emax2],..] | |
| 557 | + Keywords: | |
| 558 | + """ | |
| 559 | + #Fix all parameters except prefactor | |
| 560 | + for model in obs.models(): | |
| 561 | + for par in range(model.size()): | |
| 562 | + model[par].fix() | |
| 563 | + | |
| 564 | + opt_spectral = obs.models()[source].spectral() | |
| 565 | + opt_spectral["Prefactor"].free() | |
| 566 | + | |
| 567 | + #Define spectrum dictionary, units are in format readable by astropy | |
| 568 | + pars = {"value":[],"eu_value":[],"ed_value":[],"unit":[]} | |
| 569 | + spec = {"energy":pars.copy(),"flux":pars.copy(),"TS":[]} | |
| 570 | + | |
| 571 | + #Loop over energy bins, fit and write result to spectrum | |
| 572 | + for ebin in ebins: | |
| 573 | + #Get logarithmic bin center and flux value of all energy fit | |
| 574 | + ecenter = numpy.power(10,(numpy.log10(ebin[1])+numpy.log10(ebin[0]))/2.) | |
| 575 | + allflux = opt_spectral.eval(gammalib.GEnergy(ecenter,"TeV"),gammalib.GTime(0)) | |
| 576 | + | |
| 577 | + #Clone observations and reset energy thresholds | |
| 578 | + select = ctools.ctselect(obs) | |
| 579 | + select["ra"].real(-1.0) | |
| 580 | + select["dec"].real(-1.0) | |
| 581 | + select["rad"].real(-1.0) | |
| 582 | + select["tmin"].real(0) | |
| 583 | + select["tmax"].real(0) | |
| 584 | + select["emin"].real(ebin[0]) | |
| 585 | + select["emax"].real(ebin[1]) | |
| 586 | + select.run() | |
| 587 | + obsclone = select.obs() | |
| 588 | + | |
| 589 | + #Refit and write out to spectrum | |
| 590 | + binlike=fit(obsclone,log=False, debug=False, edisp=False) | |
| 591 | + binopt_spectral = binlike.obs().models()[source].spectral() | |
| 592 | + binflux=binopt_spectral.eval(gammalib.GEnergy(ecenter,"TeV"),gammalib.GTime(0)) | |
| 593 | + if binopt_spectral["Prefactor"].value()>0: | |
| 594 | + relerr=binopt_spectral["Prefactor"].error()/binopt_spectral["Prefactor"].value() | |
| 595 | + else: | |
| 596 | + relerr=0 | |
| 597 | + | |
| 598 | + #Store and convert to TeV from MeV | |
| 599 | + spec["energy"]["value"].append(ecenter) | |
| 600 | + spec["energy"]["ed_value"].append(ecenter-ebin[0]) | |
| 601 | + spec["energy"]["eu_value"].append(ebin[1]-ecenter) | |
| 602 | + spec["energy"]["unit"].append("1e12 eV") | |
| 603 | + | |
| 604 | + spec["flux"]["value"].append(binflux*1e6) | |
| 605 | + spec["flux"]["ed_value"].append(relerr*binflux*1e6) | |
| 606 | + spec["flux"]["eu_value"].append(relerr*binflux*1e6) | |
| 607 | + spec["flux"]["unit"].append("cm-2 1e-12 eV-1") | |
| 608 | + | |
| 609 | + spec["TS"].append(binlike.obs().models()[source].ts()) | |
| 610 | + | |
| 611 | + #Return spectrum | |
| 612 | + return spec | ... | ... |
pyext/setup.py.in
| ... | ... | @@ -121,8 +121,8 @@ setup(name='@PACKAGE_NAME@', |
| 121 | 121 | author='Juergen Knoedlseder', |
| 122 | 122 | author_email='@PACKAGE_BUGREPORT@', |
| 123 | 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 | 126 | include_dirs=include_dirs, |
| 127 | 127 | library_dirs=library_dirs, |
| 128 | 128 | libraries=libraries, | ... | ... |
scripts/Makefile.am
| ... | ... | @@ -5,8 +5,8 @@ binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py |
| 5 | 5 | |
| 6 | 6 | # Scripts to be distributed in pythondir |
| 7 | 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 | 11 | # Scripts for which symbolic links will be created in bindir |
| 12 | 12 | linkscripts = cspull cssens cstsdist cscaldb \ | ... | ... |
scripts/cspull.py
scripts/cssens.py
scripts/cstsdist.py
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 | |
| 23 | -import numpy | |
| 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" | |
| 24 | 4 | |
| 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 | |
| 519 | - | |
| 520 | -# ================================================= # | |
| 521 | -# Make spectrum for one Source in a CTA Observation # | |
| 522 | -# ================================================= # | |
| 523 | - | |
| 524 | -def getlogbins(nbins=12,emin=0.1,emax=100): | |
| 525 | - """ | |
| 526 | - Returns list with logarithmiscenergy bins in | |
| 527 | - format [[emin1,emax1],[emin2,emax2],..] | |
| 528 | - | |
| 529 | - Parameters: | |
| 530 | - nbins Number of bins | |
| 531 | - emin Minimum energy | |
| 532 | - emax Maximum energy | |
| 533 | - """ | |
| 534 | - lgemin=numpy.log10(emin) | |
| 535 | - lgemax=numpy.log10(emax) | |
| 536 | - lgbins=[] | |
| 537 | - lgedges=numpy.arange(lgemin,lgemax,(lgemax-lgemin)/nbins) | |
| 538 | - lgedges=numpy.append(lgedges,[lgemax]) | |
| 539 | - for ii in range(len(lgedges)-1): | |
| 540 | - lgbins.append([lgedges[ii],lgedges[ii+1]]) | |
| 541 | - ebins=numpy.power(10,lgbins) | |
| 542 | - return ebins | |
| 543 | - | |
| 544 | -#Please note that specpoints can make spectra for any given | |
| 545 | -#spectral model, as long as it has a "Prefactor" parameter | |
| 546 | -#which scales the funtion | |
| 547 | -def specpoints(obs,source,ebins): | |
| 548 | - """ | |
| 549 | - Returns spectral points in the given binning for | |
| 550 | - the passed observations and source. | |
| 551 | - | |
| 552 | - Parameters: | |
| 553 | - obs - Pointing direction | |
| 554 | - source - Source name for which spectrum should be done | |
| 555 | - ebins - Binning passed as a list in units of TeV as | |
| 556 | - [[emin1,emax1],[emin2,emax2],..] | |
| 557 | - Keywords: | |
| 558 | - """ | |
| 559 | - #Fix all parameters except prefactor | |
| 560 | - for model in obs.models(): | |
| 561 | - for par in range(model.size()): | |
| 562 | - model[par].fix() | |
| 563 | - | |
| 564 | - opt_spectral = obs.models()[source].spectral() | |
| 565 | - opt_spectral["Prefactor"].free() | |
| 566 | - | |
| 567 | - #Define spectrum dictionary, units are in format readable by astropy | |
| 568 | - pars = {"value":[],"eu_value":[],"ed_value":[],"unit":[]} | |
| 569 | - spec = {"energy":pars.copy(),"flux":pars.copy(),"TS":[]} | |
| 570 | - | |
| 571 | - #Loop over energy bins, fit and write result to spectrum | |
| 572 | - for ebin in ebins: | |
| 573 | - #Get logarithmic bin center and flux value of all energy fit | |
| 574 | - ecenter = numpy.power(10,(numpy.log10(ebin[1])+numpy.log10(ebin[0]))/2.) | |
| 575 | - allflux = opt_spectral.eval(gammalib.GEnergy(ecenter,"TeV"),gammalib.GTime(0)) | |
| 576 | - | |
| 577 | - #Clone observations and reset energy thresholds | |
| 578 | - select = ctools.ctselect(obs) | |
| 579 | - select["ra"].real(-1.0) | |
| 580 | - select["dec"].real(-1.0) | |
| 581 | - select["rad"].real(-1.0) | |
| 582 | - select["tmin"].real(0) | |
| 583 | - select["tmax"].real(0) | |
| 584 | - select["emin"].real(ebin[0]) | |
| 585 | - select["emax"].real(ebin[1]) | |
| 586 | - select.run() | |
| 587 | - obsclone = select.obs() | |
| 588 | - | |
| 589 | - #Refit and write out to spectrum | |
| 590 | - binlike=fit(obsclone,log=False, debug=False, edisp=False) | |
| 591 | - binopt_spectral = binlike.obs().models()[source].spectral() | |
| 592 | - binflux=binopt_spectral.eval(gammalib.GEnergy(ecenter,"TeV"),gammalib.GTime(0)) | |
| 593 | - if binopt_spectral["Prefactor"].value()>0: | |
| 594 | - relerr=binopt_spectral["Prefactor"].error()/binopt_spectral["Prefactor"].value() | |
| 595 | - else: | |
| 596 | - relerr=0 | |
| 597 | - | |
| 598 | - #Store and convert to TeV from MeV | |
| 599 | - spec["energy"]["value"].append(ecenter) | |
| 600 | - spec["energy"]["ed_value"].append(ecenter-ebin[0]) | |
| 601 | - spec["energy"]["eu_value"].append(ebin[1]-ecenter) | |
| 602 | - spec["energy"]["unit"].append("1e12 eV") | |
| 603 | - | |
| 604 | - spec["flux"]["value"].append(binflux*1e6) | |
| 605 | - spec["flux"]["ed_value"].append(relerr*binflux*1e6) | |
| 606 | - spec["flux"]["eu_value"].append(relerr*binflux*1e6) | |
| 607 | - spec["flux"]["unit"].append("cm-2 1e-12 eV-1") | |
| 608 | - | |
| 609 | - spec["TS"].append(binlike.obs().models()[source].ts()) | |
| 610 | - | |
| 611 | - #Return spectrum | |
| 612 | - return spec | |
| 5 | +from ctools.obsutils import * | ... | ... |
src/Makefile.am
| ... | ... | @@ -10,7 +10,8 @@ SUBDIRS = ctbin \ |
| 10 | 10 | ctmodel \ |
| 11 | 11 | ctobssim \ |
| 12 | 12 | ctselect \ |
| 13 | - ctskymap | |
| 13 | + ctskymap \ | |
| 14 | + ctexpcube | |
| 14 | 15 | |
| 15 | 16 | # Scripts to be distributed |
| 16 | 17 | scripts = ctools-setup \ |
| ... | ... | @@ -27,7 +28,8 @@ libctools_la_LIBADD = ctbin/libctbin.la \ |
| 27 | 28 | ctmodel/libctmodel.la \ |
| 28 | 29 | ctobssim/libctobssim.la \ |
| 29 | 30 | ctselect/libctselect.la \ |
| 30 | - ctskymap/libctskymap.la | |
| 31 | + ctskymap/libctskymap.la \ | |
| 32 | + ctexpcube/libctexpcube.la | |
| 31 | 33 | |
| 32 | 34 | # Distribute scripts |
| 33 | 35 | EXTRA_DIST = $(scripts) | ... | ... |
src/ctbin/ctbin.cpp
| ... | ... | @@ -358,15 +358,19 @@ void ctbin::save(void) |
| 358 | 358 | } |
| 359 | 359 | } |
| 360 | 360 | |
| 361 | + /* | |
| 361 | 362 | // Case A: Save counts map(s) and XML metadata information |
| 362 | 363 | if (m_use_xml) { |
| 363 | 364 | save_xml(); |
| 364 | 365 | } |
| 365 | - | |
| 366 | 366 | // Case B: Save counts map as FITS file |
| 367 | 367 | else { |
| 368 | 368 | save_fits(); |
| 369 | 369 | } |
| 370 | + */ | |
| 371 | + | |
| 372 | + // Always save FITS and an XML summary | |
| 373 | + save_xml(); | |
| 370 | 374 | |
| 371 | 375 | // Return |
| 372 | 376 | return; |
| ... | ... | @@ -757,21 +761,32 @@ void ctbin::free_members(void) |
| 757 | 761 | /***********************************************************************//** |
| 758 | 762 | * @brief Set output file name. |
| 759 | 763 | * |
| 760 | - * @param[in] filename Input file name. | |
| 764 | + * @param[in] index Index of input file. | |
| 761 | 765 | * |
| 762 | 766 | * Converts an input filename into an output filename by prepending the |
| 763 | 767 | * prefix stored in the member m_prefix to the input filename. Any path will |
| 764 | 768 | * be stripped from the input filename. Also a trailing ".gz" will be |
| 765 | - * stripped. | |
| 769 | + * stripped. If no input file name exist, prefix and observation | |
| 770 | + * number are used. | |
| 766 | 771 | ***************************************************************************/ |
| 767 | -std::string ctbin::set_outfile_name(const std::string& filename) const | |
| 772 | +std::string ctbin::set_outfile_name(const int index) const | |
| 768 | 773 | { |
| 774 | + | |
| 775 | + // Initialize output file name | |
| 776 | + std::string outname=""; | |
| 777 | + | |
| 769 | 778 | // Split input filename into path elements |
| 770 | - std::vector<std::string> elements = gammalib::split(filename, "/"); | |
| 771 | - | |
| 772 | - // The last path element is the filename | |
| 773 | - std::string outname = m_prefix + elements[elements.size()-1]; | |
| 774 | - | |
| 779 | + std::vector<std::string> elements = gammalib::split(m_infiles[index], "/"); | |
| 780 | + | |
| 781 | + // If there is an observation input file name | |
| 782 | + if (elements.size() > 0) { | |
| 783 | + outname = m_prefix + elements[elements.size()-1]; | |
| 784 | + } | |
| 785 | + // ...else, use observation number | |
| 786 | + else { | |
| 787 | + outname = m_prefix + m_obs[index]->id() + ".fits"; | |
| 788 | + } | |
| 789 | + | |
| 775 | 790 | // Strip any ".gz" |
| 776 | 791 | outname = gammalib::strip_chars(outname, ".gz"); |
| 777 | 792 | |
| ... | ... | @@ -841,7 +856,7 @@ void ctbin::save_xml(void) |
| 841 | 856 | if (obs != NULL) { |
| 842 | 857 | |
| 843 | 858 | // Set event output file name |
| 844 | - std::string outfile = set_outfile_name(m_infiles[i]); | |
| 859 | + std::string outfile = set_outfile_name(i); | |
| 845 | 860 | |
| 846 | 861 | // Store output file name in observation |
| 847 | 862 | obs->eventfile(outfile); | ... | ... |
src/ctbin/ctbin.hpp
| ... | ... | @@ -73,7 +73,7 @@ protected: |
| 73 | 73 | void init_members(void); |
| 74 | 74 | void copy_members(const ctbin& app); |
| 75 | 75 | void free_members(void); |
| 76 | - std::string set_outfile_name(const std::string& filename) const; | |
| 76 | + std::string set_outfile_name(const int index) const; | |
| 77 | 77 | void save_fits(void); |
| 78 | 78 | void save_xml(void); |
| 79 | 79 | void save_counts_map(const GCTAObservation* obs, | ... | ... |
src/ctexpcube/Makefile.am
0 → 100644
| 1 | +## Process this file with automake to produce Makefile.in | |
| 2 | + | |
| 3 | +# Program name | |
| 4 | +bin_PROGRAMS = ctexpcube | |
| 5 | + | |
| 6 | +# Sources | |
| 7 | +sources = ctexpcube.cpp | |
| 8 | + | |
| 9 | +# Libtool library | |
| 10 | +noinst_LTLIBRARIES = libctexpcube.la | |
| 11 | +libctexpcube_la_SOURCES = $(sources) | |
| 12 | + | |
| 13 | +# Program definition | |
| 14 | +ctexpcube_SOURCES = main.cpp | |
| 15 | +ctexpcube_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@ | |
| 16 | +ctexpcube_CPPFLAGS = @CPPFLAGS@ | |
| 17 | +ctexpcube_LDADD = libctexpcube.la | |
| 18 | + | |
| 19 | +# Parameter files | |
| 20 | +pfilesdir = $(prefix)/syspfiles | |
| 21 | +pfiles_DATA = ctexpcube.par | ... | ... |
src/ctexpcube/ctexpcube.cpp
0 → 100644
| 1 | +/*************************************************************************** | |
| 2 | + * ctexpcube - CTA exposure cube tool * | |
| 3 | + * ----------------------------------------------------------------------- * | |
| 4 | + * copyright (C) 2014 by Juergen Knoedlseder * | |
| 5 | + * ----------------------------------------------------------------------- * | |
| 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 | +/** | |
| 22 | + * @file ctexpcube.cpp | |
| 23 | + * @brief CTA exposure cube tool implementation | |
| 24 | + * @author Juergen Knoedlseder | |
| 25 | + */ | |
| 26 | + | |
| 27 | +/* __ Includes ___________________________________________________________ */ | |
| 28 | +#ifdef HAVE_CONFIG_H | |
| 29 | +#include <config.h> | |
| 30 | +#endif | |
| 31 | +#include <cstdio> | |
| 32 | +#include "ctexpcube.hpp" | |
| 33 | +#include "GTools.hpp" | |
| 34 | +#include "GWcs.hpp" | |
| 35 | + | |
| 36 | +/* __ Method name definitions ____________________________________________ */ | |
| 37 | +#define G_GET_EBOUNDS "ctexpcube::get_ebounds()" | |
| 38 | +#define G_GET_PARAMETERS "ctexpcube::get_parameters()" | |
| 39 | +#define G_SET_FROM_CNTMAP "ctexpcube::set_from_cntmap(std::string&)" | |
| 40 | + | |
| 41 | +/* __ Debug definitions __________________________________________________ */ | |
| 42 | + | |
| 43 | +/* __ Coding definitions _________________________________________________ */ | |
| 44 | + | |
| 45 | + | |
| 46 | +/*========================================================================== | |
| 47 | + = = | |
| 48 | + = Constructors/destructors = | |
| 49 | + = = | |
| 50 | + ==========================================================================*/ | |
| 51 | + | |
| 52 | +/***********************************************************************//** | |
| 53 | + * @brief Void constructor | |
| 54 | + ***************************************************************************/ | |
| 55 | +ctexpcube::ctexpcube(void) : GApplication(CTEXPCUBE_NAME, CTEXPCUBE_VERSION) | |
| 56 | +{ | |
| 57 | + // Initialise members | |
| 58 | + init_members(); | |
| 59 | + | |
| 60 | + // Write header into logger | |
| 61 | + log_header(); | |
| 62 | + | |
| 63 | + // Return | |
| 64 | + return; | |
| 65 | +} | |
| 66 | + | |
| 67 | + | |
| 68 | +/***********************************************************************//** | |
| 69 | + * @brief Observations constructor | |
| 70 | + * | |
| 71 | + * This method creates an instance of the class by copying an existing | |
| 72 | + * observations container. | |
| 73 | + ***************************************************************************/ | |
| 74 | +ctexpcube::ctexpcube(GObservations obs) : GApplication(CTEXPCUBE_NAME, CTEXPCUBE_VERSION) | |
| 75 | +{ | |
| 76 | + // Initialise members | |
| 77 | + init_members(); | |
| 78 | + | |
| 79 | + // Set observations | |
| 80 | + m_obs = obs; | |
| 81 | + | |
| 82 | + // Write header into logger | |
| 83 | + log_header(); | |
| 84 | + | |
| 85 | + // Return | |
| 86 | + return; | |
| 87 | +} | |
| 88 | + | |
| 89 | + | |
| 90 | + | |
| 91 | +/***********************************************************************//** | |
| 92 | + * @brief Command line constructor | |
| 93 | + * | |
| 94 | + * @param[in] argc Number of arguments in command line. | |
| 95 | + * @param[in] argv Array of command line arguments. | |
| 96 | + ***************************************************************************/ | |
| 97 | +ctexpcube::ctexpcube(int argc, char *argv[]) : | |
| 98 | + GApplication(CTEXPCUBE_NAME, CTEXPCUBE_VERSION, argc, argv) | |
| 99 | +{ | |
| 100 | + // Initialise members | |
| 101 | + init_members(); | |
| 102 | + | |
| 103 | + // Write header into logger | |
| 104 | + log_header(); | |
| 105 | + | |
| 106 | + // Return | |
| 107 | + return; | |
| 108 | +} | |
| 109 | + | |
| 110 | + | |
| 111 | +/***********************************************************************//** | |
| 112 | + * @brief Copy constructor | |
| 113 | + * | |
| 114 | + * @param[in] app ctexpcube application. | |
| 115 | + ***************************************************************************/ | |
| 116 | +ctexpcube::ctexpcube(const ctexpcube& app) : GApplication(app) | |
| 117 | +{ | |
| 118 | + // Initialise members | |
| 119 | + init_members(); | |
| 120 | + | |
| 121 | + // Copy members | |
| 122 | + copy_members(app); | |
| 123 | + | |
| 124 | + // Return | |
| 125 | + return; | |
| 126 | +} | |
| 127 | + | |
| 128 | + | |
| 129 | +/***********************************************************************//** | |
| 130 | + * @brief Destructor | |
| 131 | + ***************************************************************************/ | |
| 132 | +ctexpcube::~ctexpcube(void) | |
| 133 | +{ | |
| 134 | + // Free members | |
| 135 | + free_members(); | |
| 136 | + | |
| 137 | + // Return | |
| 138 | + return; | |
| 139 | +} | |
| 140 | + | |
| 141 | + | |
| 142 | +/*========================================================================== | |
| 143 | + = = | |
| 144 | + = Operators = | |
| 145 | + = = | |
| 146 | + ==========================================================================*/ | |
| 147 | + | |
| 148 | +/***********************************************************************//** | |
| 149 | + * @brief Assignment operator | |
| 150 | + * | |
| 151 | + * @param[in] app ctexpcube application. | |
| 152 | + * @return Returns ctexpcube application. | |
| 153 | + ***************************************************************************/ | |
| 154 | +ctexpcube& ctexpcube::operator=(const ctexpcube& app) | |
| 155 | +{ | |
| 156 | + // Execute only if object is not identical | |
| 157 | + if (this != &app) { | |
| 158 | + | |
| 159 | + // Copy base class members | |
| 160 | + this->GApplication::operator=(app); | |
| 161 | + | |
| 162 | + // Free members | |
| 163 | + free_members(); | |
| 164 | + | |
| 165 | + // Initialise members | |
| 166 | + init_members(); | |
| 167 | + | |
| 168 | + // Copy members | |
| 169 | + copy_members(app); | |
| 170 | + | |
| 171 | + } // endif: object was not identical | |
| 172 | + | |
| 173 | + // Return this object | |
| 174 | + return *this; | |
| 175 | +} | |
| 176 | + | |
| 177 | + | |
| 178 | +/*========================================================================== | |
| 179 | + = = | |
| 180 | + = Public methods = | |
| 181 | + = = | |
| 182 | + ==========================================================================*/ | |
| 183 | + | |
| 184 | +/***********************************************************************//** | |
| 185 | + * @brief Clear instance | |
| 186 | + ***************************************************************************/ | |
| 187 | +void ctexpcube::clear(void) | |
| 188 | +{ | |
| 189 | + // Free members | |
| 190 | + free_members(); | |
| 191 | + this->GApplication::free_members(); | |
| 192 | + | |
| 193 | + // Initialise members | |
| 194 | + this->GApplication::init_members(); | |
| 195 | + init_members(); | |
| 196 | + | |
| 197 | + // Return | |
| 198 | + return; | |
| 199 | +} | |
| 200 | + | |
| 201 | + | |
| 202 | +/***********************************************************************//** | |
| 203 | + * @brief Execute application | |
| 204 | + * | |
| 205 | + * This is the main execution method of the ctexpcube class. It is invoked | |
| 206 | + * when the executable is called from command line. The method generates | |
| 207 | + * the exposure cube and saves the result. | |
| 208 | + ***************************************************************************/ | |
| 209 | +void ctexpcube::execute(void) | |
| 210 | +{ | |
| 211 | + // Signal that some parameters should be read ahead | |
| 212 | + m_read_ahead = true; | |
| 213 | + | |
| 214 | + // Create the exposure cube | |
| 215 | + run(); | |
| 216 | + | |
| 217 | + // Save the exposure cube into the output FITS file | |
| 218 | + save(); | |
| 219 | + | |
| 220 | + // Return | |
| 221 | + return; | |
| 222 | +} | |
| 223 | + | |
| 224 | + | |
| 225 | +/***********************************************************************//** | |
| 226 | + * @brief Generate the model map(s) | |
| 227 | + * | |
| 228 | + * This method reads the task parameters from the parfile, sets up the | |
| 229 | + * observation container, loops over all CTA observations in the container | |
| 230 | + * and generates an exposure cube from the CTA observations. | |
| 231 | + ***************************************************************************/ | |
| 232 | +void ctexpcube::run(void) | |
| 233 | +{ | |
| 234 | + // If we're in debug mode then all output is also dumped on the screen | |
| 235 | + if (logDebug()) { | |
| 236 | + log.cout(true); | |
| 237 | + } | |
| 238 | + | |
| 239 | + // Get task parameters | |
| 240 | + get_parameters(); | |
| 241 | + | |
| 242 | + // Write parameters into logger | |
| 243 | + if (logTerse()) { | |
| 244 | + log_parameters(); | |
| 245 | + log << std::endl; | |
| 246 | + } | |
| 247 | + | |
| 248 | + // Set energy dispersion flag for all CTA observations and save old | |
| 249 | + // values in save_edisp vector | |
| 250 | + std::vector<bool> save_edisp; | |
| 251 | + save_edisp.assign(m_obs.size(), false); | |
| 252 | + for (int i = 0; i < m_obs.size(); ++i) { | |
| 253 | + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]); | |
| 254 | + if (obs != NULL) { | |
| 255 | + save_edisp[i] = obs->response().apply_edisp(); | |
| 256 | + obs->response().apply_edisp(m_apply_edisp); | |
| 257 | + } | |
| 258 | + } | |
| 259 | + | |
| 260 | + // Write observation(s) into logger | |
| 261 | + if (logTerse()) { | |
| 262 | + log << std::endl; | |
| 263 | + if (m_obs.size() > 1) { | |
| 264 | + log.header1("Observations"); | |
| 265 | + } | |
| 266 | + else { | |
| 267 | + log.header1("Observation"); | |
| 268 | + } | |
| 269 | + log << m_obs << std::endl; | |
| 270 | + } | |
| 271 | + | |
| 272 | + // Write header | |
| 273 | + if (logTerse()) { | |
| 274 | + log << std::endl; | |
| 275 | + log.header1("Generate exposure cube"); | |
| 276 | + } | |
| 277 | + | |
| 278 | + // Fill exposure | |
| 279 | + m_expcube.fill(m_obs); | |
| 280 | + | |
| 281 | + // Restore energy dispersion flag for all CTA observations | |
| 282 | + for (int i = 0; i < m_obs.size(); ++i) { | |
| 283 | + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]); | |
| 284 | + if (obs != NULL) { | |
| 285 | + obs->response().apply_edisp(save_edisp[i]); | |
| 286 | + } | |
| 287 | + } | |
| 288 | + | |
| 289 | + // Return | |
| 290 | + return; | |
| 291 | +} | |
| 292 | + | |
| 293 | + | |
| 294 | +/***********************************************************************//** | |
| 295 | + * @brief Save exposure cube | |
| 296 | + * | |
| 297 | + ***************************************************************************/ | |
| 298 | +void ctexpcube::save(void) | |
| 299 | +{ | |
| 300 | + // Write header | |
| 301 | + if (logTerse()) { | |
| 302 | + log << std::endl; | |
| 303 | + log.header1("Save exposure cube"); | |
| 304 | + } | |
| 305 | + | |
| 306 | + // Get output filename | |
| 307 | + m_outfile = (*this)["outfile"].filename(); | |
| 308 | + | |
| 309 | + // Save exposure cube | |
| 310 | + m_expcube.save(m_outfile, clobber()); | |
| 311 | + | |
| 312 | + // Return | |
| 313 | + return; | |
| 314 | +} | |
| 315 | + | |
| 316 | + | |
| 317 | +/***********************************************************************//** | |
| 318 | + * @brief Get application parameters | |
| 319 | + * | |
| 320 | + * Get all task parameters from parameter file or (if required) by querying | |
| 321 | + * the user. The parameters are read in the correct order. | |
| 322 | + * | |
| 323 | + * @todo Setup exposure cube from counts map | |
| 324 | + ***************************************************************************/ | |
| 325 | +void ctexpcube::get_parameters(void) | |
| 326 | +{ | |
| 327 | + // If we do not have any observations in the container then get an | |
| 328 | + // input file name or observation descriptor file | |
| 329 | + if (m_obs.size() == 0) { | |
| 330 | + get_obs(); | |
| 331 | + } | |
| 332 | + | |
| 333 | + // Make sure that response is set | |
| 334 | + set_response(); | |
| 335 | + | |
| 336 | + // Read energy dispersion flag | |
| 337 | + m_apply_edisp = (*this)["edisp"].boolean(); | |
| 338 | + | |
| 339 | + // If no counts map is specified then setup the exposure cube from | |
| 340 | + // the user parameters | |
| 341 | + std::string cntmap = (*this)["cntmap"].filename(); | |
| 342 | + if ((cntmap == "NONE") || (gammalib::strip_whitespace(cntmap) == "")) { | |
| 343 | + | |
| 344 | + // Get user parameters for counts map definition | |
| 345 | + std::string wcs = (*this)["proj"].string(); | |
| 346 | + std::string coordsys = (*this)["coordsys"].string(); | |
| 347 | + double xref = (*this)["xref"].real(); | |
| 348 | + double yref = (*this)["yref"].real(); | |
| 349 | + double binsz = (*this)["binsz"].real(); | |
| 350 | + int nxpix = (*this)["nxpix"].integer(); | |
| 351 | + int nypix = (*this)["nypix"].integer(); | |
| 352 | + | |
| 353 | + // Get energy definition | |
| 354 | + get_ebounds(); | |
| 355 | + | |
| 356 | + // Define exposure cube | |
| 357 | + m_expcube = GCTAExposure(wcs, coordsys, xref, yref, | |
| 358 | + -binsz, binsz, nxpix, nypix, | |
| 359 | + m_ebounds); | |
| 360 | + } | |
| 361 | + | |
| 362 | + // ... otherwise setup the exposure cube from the counts map | |
| 363 | + else { | |
| 364 | + | |
| 365 | + // Set exposure cube from counts map | |
| 366 | + set_from_cntmap(cntmap); | |
| 367 | + | |
| 368 | + } | |
| 369 | + | |
| 370 | + // Read output filename (if needed) | |
| 371 | + if (m_read_ahead) { | |
| 372 | + m_outfile = (*this)["outfile"].filename(); | |
| 373 | + } | |
| 374 | + | |
| 375 | + // Return | |
| 376 | + return; | |
| 377 | +} | |
| 378 | + | |
| 379 | + | |
| 380 | +/***********************************************************************//** | |
| 381 | + * @brief Get observation definition | |
| 382 | + * | |
| 383 | + * Get observation definition from the user parameters. | |
| 384 | + ***************************************************************************/ | |
| 385 | +void ctexpcube::get_obs(void) | |
| 386 | +{ | |
| 387 | + // Get input filename | |
| 388 | + std::string filename = (*this)["infile"].filename(); | |
| 389 | + | |
| 390 | + // Try first to open as FITS file | |
| 391 | + try { | |
| 392 | + | |
| 393 | + // Allocate CTA observation | |
| 394 | + GCTAObservation obs; | |
| 395 | + | |
| 396 | + // Load input file in CTA observation | |
| 397 | + obs.load(filename); | |
| 398 | + | |
| 399 | + // Append CTA observation to container | |
| 400 | + m_obs.append(obs); | |
| 401 | + | |
| 402 | + | |
| 403 | + } | |
| 404 | + | |
| 405 | + // ... otherwise try to open as XML file | |
| 406 | + catch (GException::fits_open_error &e) { | |
| 407 | + | |
| 408 | + // Load observations from XML file | |
| 409 | + m_obs.load(filename); | |
| 410 | + | |
| 411 | + } | |
| 412 | + | |
| 413 | + // Return | |
| 414 | + return; | |
| 415 | +} | |
| 416 | + | |
| 417 | + | |
| 418 | +/***********************************************************************//** | |
| 419 | + * @brief Set observation response | |
| 420 | + * | |
| 421 | + * Set response for all observations that so have no response. | |
| 422 | + ***************************************************************************/ | |
| 423 | +void ctexpcube::set_response(void) | |
| 424 | +{ | |
| 425 | + // Loop over all observations | |
| 426 | + for (int i = 0; i < m_obs.size(); ++i) { | |
| 427 | + | |
| 428 | + // Is this observation a CTA observation? | |
| 429 | + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]); | |
| 430 | + | |
| 431 | + // Yes ... | |
| 432 | + if (obs != NULL) { | |
| 433 | + | |
| 434 | + // Set response if it isn't set already | |
| 435 | + if (obs->response().aeff() == NULL) { | |
| 436 | + | |
| 437 | + // Load response information | |
| 438 | + std::string database = (*this)["caldb"].string(); | |
| 439 | + std::string irf = (*this)["irf"].string(); | |
| 440 | + | |
| 441 | + // Set calibration database. If specified parameter is a | |
| 442 | + // directory then use this as the pathname to the calibration | |
| 443 | + // database. Otherwise interpret this as the instrument name, | |
| 444 | + // the mission being "cta" | |
| 445 | + GCaldb caldb; | |
| 446 | + if (gammalib::dir_exists(database)) { | |
| 447 | + caldb.rootdir(database); | |
| 448 | + } | |
| 449 | + else { | |
| 450 | + caldb.open("cta", database); | |
| 451 | + } | |
| 452 | + | |
| 453 | + // Set reponse | |
| 454 | + obs->response(irf, caldb); | |
| 455 | + | |
| 456 | + } // endif: observation already has a response | |
| 457 | + | |
| 458 | + } // endif: observation was a CTA observation | |
| 459 | + | |
| 460 | + } // endfor: looped over all observations | |
| 461 | + | |
| 462 | + // Return | |
| 463 | + return; | |
| 464 | +} | |
| 465 | + | |
| 466 | + | |
| 467 | +/***********************************************************************//** | |
| 468 | + * @brief Get energy boundaries from parameters | |
| 469 | + * | |
| 470 | + * @exception GException::invalid_value | |
| 471 | + * Invalid extension name encountered. | |
| 472 | + * | |
| 473 | + * Get the energy boundaries from the user parameters. | |
| 474 | + ***************************************************************************/ | |
| 475 | +void ctexpcube::get_ebounds(void) | |
| 476 | +{ | |
| 477 | + // Determine the energy binning alogrithm | |
| 478 | + std::string ebinalg = (*this)["ebinalg"].string(); | |
| 479 | + | |
| 480 | + // If we have the binning given by a file then try to get the boundaries | |
| 481 | + // from that file | |
| 482 | + if (ebinalg == "FILE") { | |
| 483 | + | |
| 484 | + // Get filename | |
| 485 | + std::string filename = (*this)["ebinfile"].filename(); | |
| 486 | + | |
| 487 | + // Open fits file to check which extension is given | |
| 488 | + GFits file(filename); | |
| 489 | + | |
| 490 | + // Check first for EBOUNDS extension | |
| 491 | + if (file.contains("EBOUNDS")) { | |
| 492 | + file.close(); | |
| 493 | + m_ebounds.load(filename, "EBOUNDS"); | |
| 494 | + } | |
| 495 | + | |
| 496 | + // ... then check for ENERGYBINS extension | |
| 497 | + else if (file.contains("ENERGYBINS")) { | |
| 498 | + file.close(); | |
| 499 | + m_ebounds.load(filename, "ENERGYBINS"); | |
| 500 | + } | |
| 501 | + | |
| 502 | + // ... otherwise throw an exception | |
| 503 | + else { | |
| 504 | + file.close(); | |
| 505 | + std::string msg = "No extension with name \"EBOUNDS\" or" | |
| 506 | + " \"ENERGYBINS\" found in FITS file" | |
| 507 | + " \""+filename+"\".\n" | |
| 508 | + "An \"EBOUNDS\" or \"ENERGYBINS\" extension" | |
| 509 | + " is required if the parameter \"ebinalg\"" | |
| 510 | + " is set to \"FILE\"."; | |
| 511 | + throw GException::invalid_value(G_GET_EBOUNDS, msg); | |
| 512 | + } | |
| 513 | + } | |
| 514 | + | |
| 515 | + // ... otherwise read emin, emax and nebins | |
| 516 | + else { | |
| 517 | + | |
| 518 | + // Get the relevant parameters | |
| 519 | + double emin = (*this)["emin"].real(); | |
| 520 | + double emax = (*this)["emax"].real(); | |
| 521 | + int enumbins = (*this)["enumbins"].integer(); | |
| 522 | + bool log = ((*this)["ebinalg"].string() == "LIN") ? false : true; | |
| 523 | + | |
| 524 | + // Create energy boundaries | |
| 525 | + m_ebounds = GEbounds(enumbins, | |
| 526 | + GEnergy(emin, "TeV"), | |
| 527 | + GEnergy(emax, "TeV"), | |
| 528 | + log); | |
| 529 | + | |
| 530 | + } | |
| 531 | + | |
| 532 | + // Return | |
| 533 | + return; | |
| 534 | +} | |
| 535 | + | |
| 536 | + | |
| 537 | +/***********************************************************************//** | |
| 538 | + * @brief Set exposure cube definition from counts map | |
| 539 | + * | |
| 540 | + * @exception GException::invalid_value | |
| 541 | + * Invalid counts map projection or invalid events encountered. | |
| 542 | + * | |
| 543 | + * Set exposure cube definition from counts map. | |
| 544 | + ***************************************************************************/ | |
| 545 | +void ctexpcube::set_from_cntmap(const std::string& filename) | |
| 546 | +{ | |
| 547 | + // Allocate CTA observation | |
| 548 | + GCTAObservation obs; | |
| 549 | + | |
| 550 | + // Load counts map in CTA observation | |
| 551 | + obs.load(filename); | |
| 552 | + | |
| 553 | + // Set exposure cube from counts map | |
| 554 | + const GCTAEventCube* cube = dynamic_cast<const GCTAEventCube*>(obs.events()); | |
| 555 | + | |
| 556 | + // Continue only if cube is valid | |
| 557 | + if (cube != NULL) { | |
| 558 | + | |
| 559 | + // Get sky map projection | |
| 560 | + const GWcs* wcs = dynamic_cast<const GWcs*>(cube->map().projection()); | |
| 561 | + | |
| 562 | + // Continue only if projection is valid | |
| 563 | + if (wcs != NULL) { | |
| 564 | + | |
| 565 | + // Get user parameters for counts map definition | |
| 566 | + std::string proj = wcs->code(); | |
| 567 | + std::string coordsys = wcs->coordsys(); | |
| 568 | + double xref = wcs->crval(0); | |
| 569 | + double yref = wcs->crval(1); | |
| 570 | + double dx = wcs->cdelt(0); | |
| 571 | + double dy = wcs->cdelt(1); | |
| 572 | + int nx = cube->map().nx(); | |
| 573 | + int ny = cube->map().ny(); | |
| 574 | + | |
| 575 | + // Get energy definition | |
| 576 | + m_ebounds = cube->ebounds(); | |
| 577 | + | |
| 578 | + // Define exposure cube | |
| 579 | + m_expcube = GCTAExposure(proj, coordsys, xref, yref, | |
| 580 | + dx, dy, nx, ny, | |
| 581 | + m_ebounds); | |
| 582 | + | |
| 583 | + } // endif: WCS projection was valid | |
| 584 | + | |
| 585 | + // ... projection is not of WCS type | |
| 586 | + else { | |
| 587 | + std::string msg = "Counts map project is not of WCS type."; | |
| 588 | + throw GException::invalid_value(G_SET_FROM_CNTMAP, msg); | |
| 589 | + } | |
| 590 | + | |
| 591 | + } // endif: observation contained an events cube | |
| 592 | + | |
| 593 | + // ... there is not events cube | |
| 594 | + else { | |
| 595 | + std::string msg = "No events cube found in file \"" | |
| 596 | + ""+filename+"\"."; | |
| 597 | + throw GException::invalid_value(G_SET_FROM_CNTMAP, msg); | |
| 598 | + } | |
| 599 | + | |
| 600 | + // Return | |
| 601 | + return; | |
| 602 | +} | |
| 603 | + | |
| 604 | + | |
| 605 | +/*========================================================================== | |
| 606 | + = = | |
| 607 | + = Private methods = | |
| 608 | + = = | |
| 609 | + ==========================================================================*/ | |
| 610 | + | |
| 611 | +/***********************************************************************//** | |
| 612 | + * @brief Initialise class members | |
| 613 | + ***************************************************************************/ | |
| 614 | +void ctexpcube::init_members(void) | |
| 615 | +{ | |
| 616 | + // Initialise members | |
| 617 | + m_outfile.clear(); | |
| 618 | + m_apply_edisp = false; | |
| 619 | + | |
| 620 | + // Initialise protected members | |
| 621 | + m_read_ahead = false; | |
| 622 | + m_obs.clear(); | |
| 623 | + m_expcube.clear(); | |
| 624 | + m_ebounds.clear(); | |
| 625 | + | |
| 626 | + // Set logger properties | |
| 627 | + log.date(true); | |
| 628 | + | |
| 629 | + // Return | |
| 630 | + return; | |
| 631 | +} | |
| 632 | + | |
| 633 | + | |
| 634 | +/***********************************************************************//** | |
| 635 | + * @brief Copy class members | |
| 636 | + * | |
| 637 | + * @param[in] app Application. | |
| 638 | + ***************************************************************************/ | |
| 639 | +void ctexpcube::copy_members(const ctexpcube& app) | |
| 640 | +{ | |
| 641 | + // Copy attributes | |
| 642 | + m_outfile = app.m_outfile; | |
| 643 | + m_apply_edisp = app.m_apply_edisp; | |
| 644 | + | |
| 645 | + // Copy protected members | |
| 646 | + m_read_ahead = app.m_read_ahead; | |
| 647 | + m_obs = app.m_obs; | |
| 648 | + m_expcube = app.m_expcube; | |
| 649 | + m_ebounds = app.m_ebounds; | |
| 650 | + | |
| 651 | + // Return | |
| 652 | + return; | |
| 653 | +} | |
| 654 | + | |
| 655 | + | |
| 656 | +/***********************************************************************//** | |
| 657 | + * @brief Delete class members | |
| 658 | + ***************************************************************************/ | |
| 659 | +void ctexpcube::free_members(void) | |
| 660 | +{ | |
| 661 | + // Write separator into logger | |
| 662 | + if (logTerse()) { | |
| 663 | + log << std::endl; | |
| 664 | + } | |
| 665 | + | |
| 666 | + // Return | |
| 667 | + return; | |
| 668 | +} | ... | ... |
src/ctexpcube/ctexpcube.hpp
0 → 100644
| 1 | +/*************************************************************************** | |
| 2 | + * ctexpcube - CTA exposure cube tool * | |
| 3 | + * ----------------------------------------------------------------------- * | |
| 4 | + * copyright (C) 2014 by Juergen Knoedlseder * | |
| 5 | + * ----------------------------------------------------------------------- * | |
| 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 | +/** | |
| 22 | + * @file ctexpcube.hpp | |
| 23 | + * @brief CTA exposure cube tool definition | |
| 24 | + * @author Juergen Knoedlseder | |
| 25 | + */ | |
| 26 | + | |
| 27 | +#ifndef CTEXPCUBE_HPP | |
| 28 | +#define CTEXPCUBE_HPP | |
| 29 | + | |
| 30 | +/* __ Includes ___________________________________________________________ */ | |
| 31 | +#include "GammaLib.hpp" | |
| 32 | +#include "GCTALib.hpp" | |
| 33 | + | |
| 34 | +/* __Definitions _________________________________________________________ */ | |
| 35 | +#define CTEXPCUBE_NAME "ctexpcube" | |
| 36 | +#define CTEXPCUBE_VERSION "00-01-00" | |
| 37 | + | |
| 38 | + | |
| 39 | +/***********************************************************************//** | |
| 40 | + * @class ctexpcube | |
| 41 | + * | |
| 42 | + * @brief CTA exposure cube tool | |
| 43 | + ***************************************************************************/ | |
| 44 | +class ctexpcube : public GApplication { | |
| 45 | + | |
| 46 | +public: | |
| 47 | + // Constructors and destructors | |
| 48 | + ctexpcube(void); | |
| 49 | + explicit ctexpcube(GObservations obs); | |
| 50 | + ctexpcube(int argc, char *argv[]); | |
| 51 | + ctexpcube(const ctexpcube& app); | |
| 52 | + virtual ~ctexpcube(void); | |
| 53 | + | |
| 54 | + // Operators | |
| 55 | + ctexpcube& operator=(const ctexpcube& app); | |
| 56 | + | |
| 57 | + // Methods | |
| 58 | + void clear(void); | |
| 59 | + void execute(void); | |
| 60 | + void run(void); | |
| 61 | + void save(void); | |
| 62 | + GObservations& obs(void) { return m_obs; } | |
| 63 | + void get_parameters(void); | |
| 64 | + void get_obs(void); | |
| 65 | + void set_response(void); | |
| 66 | + void get_ebounds(void); | |
| 67 | + void set_from_cntmap(const std::string& filename); | |
| 68 | + | |
| 69 | +protected: | |
| 70 | + // Protected methods | |
| 71 | + void init_members(void); | |
| 72 | + void copy_members(const ctexpcube& app); | |
| 73 | + void free_members(void); | |
| 74 | + | |
| 75 | + // User parameters | |
| 76 | + std::string m_outfile; //!< Output exposure cube file | |
| 77 | + bool m_apply_edisp; //!< Apply energy dispersion? | |
| 78 | + | |
| 79 | + // Protected members | |
| 80 | + bool m_read_ahead; //!< Read ahead parameters | |
| 81 | + GObservations m_obs; //!< Observation container | |
| 82 | + GCTAExposure m_expcube; //!< Exposure cube | |
| 83 | + GEbounds m_ebounds; //!< Energy boundaries | |
| 84 | +}; | |
| 85 | + | |
| 86 | +#endif /* CTEXPCUBE_HPP */ | ... | ... |
src/ctexpcube/ctexpcube.par
0 → 100644
| 1 | +############################################################################# | |
| 2 | +# # | |
| 3 | +# CTA Science Analysis Tools # | |
| 4 | +# # | |
| 5 | +# ------------------------------------------------------------------------- # | |
| 6 | +# # | |
| 7 | +# File: ctexpcube.par # | |
| 8 | +# Version: 00-01-00 # | |
| 9 | +# # | |
| 10 | +# Author: Juergen Knoedlseder # | |
| 11 | +# Jurgen.Knodlseder@irap.omp.eu # | |
| 12 | +# IRAP # | |
| 13 | +# # | |
| 14 | +# Purpose: Parameter file for the ctexpcube application # | |
| 15 | +# # | |
| 16 | +############################################################################# | |
| 17 | +# | |
| 18 | +# General parameters | |
| 19 | +#=================== | |
| 20 | +infile, f, a, "NONE",,, "Observation definition file" | |
| 21 | +cntmap, f, a, "NONE",,, "Counts map for exposure cube definition" | |
| 22 | +outfile, f, a, "expcube.fits",,, "Output exposure cube file" | |
| 23 | +caldb, s, a, "dummy",,, "Calibration database" | |
| 24 | +irf, s, a, "cta_dummy_irf",,, "Instrument response function" | |
| 25 | +edisp, b, h, no,,, "Apply energy dispersion?" | |
| 26 | +# | |
| 27 | +# Energy binning parameters | |
| 28 | +#========================== | |
| 29 | +ebinalg, s, h, "LOG", FILE|LIN|LOG,,"Algorithm for defining energy bins" | |
| 30 | +emin, r, a, 0.1,,,"Start value for first energy bin in TeV" | |
| 31 | +emax, r, a, 100.0,,,"Stop value for last energy bin in TeV" | |
| 32 | +enumbins, i, a, 20,,,"Number of energy bins" | |
| 33 | +ebinfile, f, h, "NONE",,,"Name of the file containing the energy bin definition" | |
| 34 | +# | |
| 35 | +# Spatial binning parameters | |
| 36 | +#=========================== | |
| 37 | +nxpix, i, a, 200,,, "Size of the X axis in pixels" | |
| 38 | +nypix, i, a, 200,,, "Size of the Y axis in pixels" | |
| 39 | +binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" | |
| 40 | +coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" | |
| 41 | +xref, r, a, 83.63,,, "First coordinate of image center in degrees (RA or galactic l)" | |
| 42 | +yref, r, a, 22.01,,, "Second coordinate of image center in degrees (DEC or galactic b)" | |
| 43 | +axisrot, r, h, 0.0,,, "Rotation angle of image axis, in degrees" | |
| 44 | +proj, s, a, "CAR", AIT|AZP|CAR|MER|STG|TAN,, "Projection method e.g. AIT|AZP|CAR|MER|STG|TAN" | |
| 45 | +# | |
| 46 | +# Standard parameters | |
| 47 | +#==================== | |
| 48 | +chatter, i, h, 2,0,4, "Chattiness of output" | |
| 49 | +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" | |
| 50 | +debug, b, h, no,,, "Debugging mode activated" | |
| 51 | +mode, s, h, "ql",,, "Mode of automatic parameters" | |
| 52 | +logfile, f, h, "ctexpcube.log",,, "Log filename" | ... | ... |
src/ctexpcube/main.cpp
0 → 100644
| 1 | +/*************************************************************************** | |
| 2 | + * ctexpcube - CTA exposure cube tool * | |
| 3 | + * ----------------------------------------------------------------------- * | |
| 4 | + * copyright (C) 2014 by Juergen Knoedlseder * | |
| 5 | + * ----------------------------------------------------------------------- * | |
| 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 | +/** | |
| 22 | + * @file ctexpcube/main.cpp | |
| 23 | + * @brief CTA exposure cube tool main code | |
| 24 | + * @author Juergen Knoedlseder | |
| 25 | + */ | |
| 26 | + | |
| 27 | +/* __ Includes ___________________________________________________________ */ | |
| 28 | +#ifdef HAVE_CONFIG_H | |
| 29 | +#include <config.h> | |
| 30 | +#endif | |
| 31 | +#include <cstdio> | |
| 32 | +#include "ctexpcube.hpp" | |
| 33 | + | |
| 34 | + | |
| 35 | +/***********************************************************************//** | |
| 36 | + * @brief Main entry point of application | |
| 37 | + * | |
| 38 | + * @param[in] argc Number of command line arguments. | |
| 39 | + * @param[in] argv Command line arguments. | |
| 40 | + * | |
| 41 | + * This is the main entry point of the ctexpcube application. It allocates a | |
| 42 | + * ctexpcube object and runs the application. | |
| 43 | + ***************************************************************************/ | |
| 44 | +int main (int argc, char *argv[]) | |
| 45 | +{ | |
| 46 | + // Initialise return code | |
| 47 | + int rc = 1; | |
| 48 | + | |
| 49 | + // Create instance of application | |
| 50 | + ctexpcube application(argc, argv); | |
| 51 | + | |
| 52 | + // Run application | |
| 53 | + try { | |
| 54 | + // Execute application | |
| 55 | + application.execute(); | |
| 56 | + | |
| 57 | + // Signal success | |
| 58 | + rc = 0; | |
| 59 | + } | |
| 60 | + catch (std::exception &e) { | |
| 61 | + | |
| 62 | + // Extract error message | |
| 63 | + std::string message = e.what(); | |
| 64 | + std::string signal = "*** ERROR encounterted in the execution of" | |
| 65 | + " ctexpcube. Run aborted ..."; | |
| 66 | + | |
| 67 | + // Write error in logger | |
| 68 | + application.log << signal << std::endl; | |
| 69 | + application.log << message << std::endl; | |
| 70 | + | |
| 71 | + // Write error on standard output | |
| 72 | + std::cout << signal << std::endl; | |
| 73 | + std::cout << message << std::endl; | |
| 74 | + | |
| 75 | + } // endcatch: catched any application error | |
| 76 | + | |
| 77 | + // Return | |
| 78 | + return rc; | |
| 79 | +} | ... | ... |
test/Makefile.am
| ... | ... | @@ -8,7 +8,7 @@ MACOSX_DEPLOYMENT_TARGET=@CONFIGURE_MACOSX_DEPLOYMENT_TARGET@ |
| 8 | 8 | if WITH_PYTHON |
| 9 | 9 | TEST_CSCRIPTS = test_cscripts.sh |
| 10 | 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 | 12 | endif |
| 13 | 13 | |
| 14 | 14 | # Set test environment | ... | ... |
test/data/crab_cntmap.fits.gz
No preview for this file type
test/data/crab_events.fits.gz
No preview for this file type
test/data/run_00006028_cntmap.fits.gz deleted
No preview for this file type
test/data/run_00006028_eventlist_reco.fits.gz deleted
No preview for this file type
test/example_survey.jpg deleted
245 KB
test/test_ctools.sh
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | # ========================================================================== |
| 3 | 3 | # This script tests all ctools |
| 4 | 4 | # |
| 5 | -# Copyright (C) 2011-2012 Juergen Knoedlseder | |
| 5 | +# Copyright (C) 2011-2014 Juergen Knoedlseder | |
| 6 | 6 | # |
| 7 | 7 | # This program is free software: you can redistribute it and/or modify |
| 8 | 8 | # it under the terms of the GNU General Public License as published by |
| ... | ... | @@ -40,6 +40,7 @@ ctmodel=../src/ctmodel/ctmodel |
| 40 | 40 | ctobssim=../src/ctobssim/ctobssim |
| 41 | 41 | ctselect=../src/ctselect/ctselect |
| 42 | 42 | ctskymap=../src/ctskymap/ctskymap |
| 43 | +ctexpcube=../src/ctexpcube/ctexpcube | |
| 43 | 44 | |
| 44 | 45 | |
| 45 | 46 | # |
| ... | ... | @@ -340,4 +341,49 @@ else |
| 340 | 341 | exit 1 |
| 341 | 342 | fi |
| 342 | 343 | $ECHO " ok" |
| 343 | -#sleep 100 | |
| 344 | + | |
| 345 | + | |
| 346 | +# | |
| 347 | +# Test ctexpcube | |
| 348 | +# ============== | |
| 349 | +$ECHO -n "Test ctexpcube: " | |
| 350 | +# | |
| 351 | +# Run 1 | |
| 352 | +$ctexpcube infile="data/crab_events.fits" \ | |
| 353 | + cntmap="NONE" \ | |
| 354 | + outfile="expcube1.fits" \ | |
| 355 | + caldb="irf" \ | |
| 356 | + irf="cta_dummy_irf" \ | |
| 357 | + emin=0.1 \ | |
| 358 | + emax=100.0 \ | |
| 359 | + enumbins=20 \ | |
| 360 | + nxpix=200 \ | |
| 361 | + nypix=200 \ | |
| 362 | + binsz=0.02 \ | |
| 363 | + coordsys="CEL" \ | |
| 364 | + xref=83.63 \ | |
| 365 | + yref=22.01 \ | |
| 366 | + proj="CAR" | |
| 367 | +$ECHO -n "." | |
| 368 | +if [ -s "expcube1.fits" ] | |
| 369 | +then | |
| 370 | + $ECHO -n "." | |
| 371 | +else | |
| 372 | + $ECHO " expcube1.fits file is not valid" | |
| 373 | + exit 1 | |
| 374 | +fi | |
| 375 | +# | |
| 376 | +# Run 2 | |
| 377 | +$ctexpcube infile="data/crab_events.fits" \ | |
| 378 | + cntmap="data/crab_cntmap.fits" \ | |
| 379 | + outfile="expcube2.fits" \ | |
| 380 | + caldb="irf" \ | |
| 381 | + irf="cta_dummy_irf" | |
| 382 | +$ECHO -n "." | |
| 383 | +if [ -s "expcube2.fits" ] | |
| 384 | +then | |
| 385 | + $ECHO -n "." | |
| 386 | +else | |
| 387 | + $ECHO " expcube2.fits file is not valid" | |
| 388 | + exit 1 | |
| 389 | +fi | ... | ... |
test/test_python.py
| ... | ... | @@ -228,7 +228,7 @@ if __name__ == '__main__': |
| 228 | 228 | |
| 229 | 229 | # Test diffuse map cube with constant spectrum |
| 230 | 230 | sys.stdout.write("Test diffuse map cube with constant analysis: ") |
| 231 | - ref = {'CTABackgroundModel': {'Value': 1.03547}} | |
| 231 | + ref = {'CTABackgroundModel': {'Value': 0.242717}} | |
| 232 | 232 | err = pipeline_v2(model_name="data/model_background_cube_const.xml", \ |
| 233 | 233 | ra_sim=84.17263, dec_sim=22.01444, rad_sim=2.0, \ |
| 234 | 234 | ra_select=83.80, dec_select=22.20, rad_select=1.0, \ |
| ... | ... | @@ -242,7 +242,7 @@ if __name__ == '__main__': |
| 242 | 242 | |
| 243 | 243 | # Test diffuse map cube with power law spectrum |
| 244 | 244 | sys.stdout.write("Test diffuse map cube with power law analysis: ") |
| 245 | - ref = {'CTABackgroundModel': {'Prefactor': 1.02396, 'Index': -0.00710087}} | |
| 245 | + ref = {'CTABackgroundModel': {'Prefactor': 0.513946, 'Index': -1.22388}} | |
| 246 | 246 | err = pipeline_v2(model_name="data/model_background_cube.xml", \ |
| 247 | 247 | ra_sim=84.17263, dec_sim=22.01444, rad_sim=2.0, \ |
| 248 | 248 | ra_select=84.00, dec_select=22.20, rad_select=1.0, \ | ... | ... |