Dmitry Malyshev / ctools

Commit 89ab90bbe67d8ce12d5553175505e87d4bdc74e7

Authored by Jürgen Knödlseder
2 parents 5f3ab77e ec7d56af

Merge branch '00-07-02' into integration 

Conflicts:
	ChangeLog
	scripts/cta_root2caldb.py
	src/ctlike/ctlike.par
	src/ctobssim/ctobssim.par
Inline Side-by-side
Showing 26 changed files with 1524 additions and 272 deletions
ChangeLog
  View file @89ab90b
... ... @@ -6,12 +6,13 @@
6 6 Support energy dispersion in ctobssim, ctlike and ctmodel
7 7  
8 8  
9   -2014-02-10 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  9 +2014-02-20 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
10 10  
11 11 * Bug fix Version 0.7.2 released
12 12 ================================
13 13  
14 14 Add logL and event statistics to cspull output file
  15 + Properly close temporary file in ctselect (#1149)
15 16  
16 17  
17 18 2014-02-05 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
... ...
... ... @@ -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   -10 February 2014
  5 +20 February 2014
6 6  
7 7  
8 8 1. Introduction
... ... @@ -31,7 +31,8 @@ dispersion can be switched on or off using the hidden edisp parameter.
31 31  
32 32 5. ctselect - CTA event selection
33 33 ---------------------------------
34   -None
  34 +The temporary file is now properly close to avoid file opening errors
  35 +when iterating over ctselect (#1149).
35 36  
36 37  
37 38 6. ctbin - CTA event binning
... ...
caldb/Makefile.am
  View file @89ab90b
1 1 ## Process this file with automake to produce Makefile.in
2 2  
3 3 # Calibration files to distribute
4   -irfs = cta_dummy_irf.dat
  4 +caldb_indx = caldb.indx
  5 +caldb_irf = cta_dummy_irf.dat
5 6  
6   -caldbdir = $(datadir)/caldb/data/cta/dummy
7   -caldb_DATA = $(irfs)
8   -
9   -# Extra files to be distributed in the tar file
10   -EXTRA_DIST = $(irfs)
  7 +# Distribute calibration database files
  8 +EXTRA_DIST = $(caldb_indx) $(caldb_irf)
  9 +caldbindxdir = $(datadir)/caldb/data/cta/dummy
  10 +caldbindx_DATA = $(caldb_indx)
  11 +caldbirfdir = $(datadir)/caldb/data/cta/dummy/bcf
  12 +caldbirf_DATA = $(caldb_irf)
... ...
caldb/caldb.indx 0 → 100644
View file @89ab90b
No preview for this file type
doc/source/_templates/index.html
  View file @89ab90b
... ... @@ -34,7 +34,8 @@ CTA format), Fermi/LAT and COMPTEL data. {%endtrans%}&lt;/p&gt;
34 34  
35 35 <ul>
36 36 <li><a class="smalllink" href="{{ pathto("user_manual/getting") }}">{%trans%}Installation{%endtrans%}</a></li>
37   -<li><a class="smalllink" href="{{ pathto("user_manual/quickstart") }}">{%trans%}Quickstart{%endtrans%}</a></li>
  37 +<li><a class="smalllink" href="{{ pathto("user_manual/getting_started/quickstart") }}">{%trans%}Quickstart{%endtrans%}</a></li>
  38 +<li><a class="smalllink" href="{{ pathto("user_manual/getting_started/response") }}">{%trans%}CTA response functions{%endtrans%}</a></li>
38 39 <li><a class="smalllink" href="{{ pathto("user_manual/index") }}">{%trans%}User Manual{%endtrans%}</a></li>
39 40 <li><a class="smalllink" href="{{ pathto("documentation") }}">{%trans%}ctools reference{%endtrans%}</a></li>
40 41 <li><a class="smalllink" href="doxygen/index.html">{%trans%}Source code documentation (Doxygen){%endtrans%}</a></li>
... ...
doc/source/conf.py
  View file @89ab90b
... ... @@ -51,7 +51,7 @@ copyright = u&#39;2014, ctools developers&#39;
51 51 version = '0.7'
52 52  
53 53 # The full version, including alpha/beta/rc tags.
54   -release = '0.7.1'
  54 +release = '0.7.2'
55 55  
56 56 # The language for content autogenerated by Sphinx. Refer to documentation
57 57 # for a list of supported languages.
... ... @@ -165,7 +165,7 @@ html_show_copyright = False
165 165 #html_file_suffix = None
166 166  
167 167 # Output file base name for HTML help builder.
168   -htmlhelp_basename = 'GammaLibdoc'
  168 +htmlhelp_basename = 'ctoolsdoc'
169 169  
170 170  
171 171 # -- Options for LaTeX output --------------------------------------------------
... ...
doc/source/user_manual/cntmap-fits.jpg renamed to doc/source/user_manual/getting_started/cntmap-fits.jpg
View file @89ab90b

47 KB
doc/source/user_manual/cntmap-map.jpg renamed to doc/source/user_manual/getting_started/cntmap-map.jpg
View file @89ab90b

22.7 KB
doc/source/user_manual/dskeys.jpg renamed to doc/source/user_manual/getting_started/dskeys.jpg
View file @89ab90b

99.6 KB
doc/source/user_manual/getting_started/faq.rst 0 → 100644
  View file @89ab90b
  1 +Frequently asked questions
  2 +--------------------------
  3 +
  4 +To be written ...
... ...
doc/source/user_manual/getting_started/index.rst 0 → 100644
  View file @89ab90b
  1 +Getting started
  2 +===============
  3 +
  4 +.. toctree::
  5 + :maxdepth: 2
  6 +
  7 + quickstart
  8 + response
  9 + models
  10 + python
  11 + faq
  12 +
... ...
doc/source/user_manual/getting_started/irfs.jpg 0 → 100644
View file @89ab90b

27.6 KB
doc/source/user_manual/getting_started/models.rst 0 → 100644
  View file @89ab90b
  1 +Models available in ctools
  2 +--------------------------
  3 +
  4 +To be written ... (you may read
  5 +http://gammalib.sourceforge.net/user_manual/modules/model.html#spectral-models for the time being).
... ...
doc/source/user_manual/getting_started/python.rst 0 → 100644
  View file @89ab90b
  1 +Using ctools from Python
  2 +------------------------
  3 +
  4 +To be written ...
... ...
doc/source/user_manual/quickstart.rst renamed to doc/source/user_manual/getting_started/quickstart.rst
  View file @89ab90b
1   -Getting started
2   -===============
3   -
4 1 A quick ctools tutorial
5 2 -----------------------
6 3  
... ... @@ -641,22 +638,3 @@ the same events for both analyses.
641 638 2014-01-09T21:47:54: +==============+
642 639 2014-01-09T21:47:54:
643 640 2014-01-09T21:47:54: Application "ctlike" terminated after 6 wall clock seconds, consuming 0.163384 seconds of CPU time.
644   -
645   -
646   -Models available in ctools
647   ---------------------------
648   -
649   -To be written ... (you may read
650   -http://gammalib.sourceforge.net/user_manual/modules/model.html#spectral-models for the time being).
651   -
652   -
653   -Using ctools from Python
654   -------------------------
655   -
656   -To be written ...
657   -
658   -
659   -Frequently asked questions
660   ---------------------------
661   -
662   -To be written ...
... ...
doc/source/user_manual/getting_started/response.rst 0 → 100644
  View file @89ab90b
  1 +CTA Instrument Response Functions
  2 +---------------------------------
  3 +
  4 +.. note ::
  5 +
  6 + Instrument response functions are formally not part of ctools as they
  7 + should be provided by the instrument teams. The GammaLib framework on
  8 + which ctools are built comes with a set of basic response functions, but
  9 + for getting the latest instrument response functions you should contact
  10 + the relevant instrument teams. CTA consortium members should
  11 + `download
  12 + <https://portal.cta-observatory.org/WG/DM/DM_wiki/DATA_Access/Pages/Science%20Tools.aspx>`_
  13 + the latest calibration database containing Prod1 and Prod2 instrument
  14 + response functions from the consortium Sharepoint site (requires
  15 + password).
  16 +
  17 +
  18 +What are instrument response functions?
  19 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  20 +
  21 +The instrument response functions provide a mathematical description that
  22 +links the measured quantities :math:`\vec{d}` of an event to the physical
  23 +quantities :math:`\vec{p}` of the incident photon. The following figure
  24 +illustrates this relationship:
  25 +
  26 +.. figure:: irfs.jpg
  27 + :width: 70%
  28 +
  29 +:math:`I(\vec{p})` is the gamma-ray intensity arriving at Earth as a
  30 +function of photon properties :math:`\vec{p}`
  31 +(which usually are true photon energy, true photon incident direction,
  32 +and true photon arrival time),
  33 +while :math:`e(\vec{d})` is the expected event rate as function of event
  34 +properties :math:`\vec{d}` (which usually are the measured photon energy,
  35 +measured or reconstructued photon incident direction, and measured photon
  36 +arrival time). The expected event rate is obtained by integrating the
  37 +product of
  38 +the instrumental response function :math:`R(\vec{d}|\vec{p},\vec{a}`)
  39 +and the emitted intensity :math:`I(\vec{p})` over the photon properties
  40 +:math:`\vec{p}`.
  41 +The argument :math:`\vec{a}` in the response function comprises any
  42 +auxiliary parameter on which the response function may depend on (e.g.
  43 +pointing direction, triggered telescopes, optical efficiencies,
  44 +atmospheric conditions, etc.). All these quantities and hence the
  45 +instrument response function may depend on time.
  46 +
  47 +
  48 +Data formats
  49 +~~~~~~~~~~~~
  50 +
  51 +So far the data format of the CTA instrument response functions has not
  52 +yet been settled, and information about instrument response functions is
  53 +provided in various ways. Three types of data formats are supported:
  54 +
  55 +- :ref:`sec_cta_rsp_perftable`
  56 +
  57 +- :ref:`sec_cta_rsp_xspec`
  58 +
  59 +- :ref:`sec_cta_rsp_rsptable`
  60 +
  61 +All formats provide the instrument response as function of true (and
  62 +sometimes also measured) photon energy, typically from about 20 GeV to
  63 +about 125 TeV.
  64 +Performance tables provide the instrument response for on-axis sources.
  65 +ARF, RMF and PSF files provide the instrument response for a specific
  66 +source position (and region) within the field of view.
  67 +Response tables provide the instrument response as function of position in
  68 +the field of view. Response tables are thus the most universal form of
  69 +instrument response functions available, and we recommend using this form
  70 +for any analysis.
  71 +
  72 +.. note ::
  73 +
  74 + Instrument response functions are so far only provided for a fixed
  75 + zenith angle of 20 deg. Therefore, no zenith or azimuth angle dependence has
  76 + been implemented so far.
  77 +
  78 +
  79 +Installing the CTA calibration database
  80 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  81 +
  82 +After `downloading
  83 +<https://portal.cta-observatory.org/WG/DM/DM_wiki/DATA_Access/Pages/Science%20Tools.aspx>`_
  84 +the latest calibration database (only possible for CTA consortium members),
  85 +the database is installed using
  86 +
  87 +.. code-block:: bash
  88 +
  89 + $ [sudo] tar -C $CALDB -zxvf cta-caldb-20140220.tar.gz
  90 +
  91 +or
  92 +
  93 +.. code-block:: bash
  94 +
  95 + $ [sudo] tar -C $CALDB -xvf cta-caldb-20140220.tar
  96 +
  97 +depending on whether the database has been retrieved as a gzipped file or
  98 +not. Depending on your access rights to the ``$CALDB`` directory,
  99 +installation of the calibration database may require root privileges
  100 +(type ``sudo`` in this case, otherwise omit this part of the command).
  101 +
  102 +The calibration database must be installed into the directory to which you
  103 +``CALDB`` environment variable points. By default, ctools sets this
  104 +enviroment variable to
  105 +
  106 +.. code-block:: bash
  107 +
  108 + $CTOOLS/share/caldb
  109 +
  110 +but it may be that due to some other analysis software on your system the
  111 +``CALDB`` environment variable points to another directory. If in doubt, type
  112 +
  113 +.. code-block:: bash
  114 +
  115 + echo $CALDB
  116 +
  117 +to find out to which directory your ``CALDB`` environment variable points.
  118 +You may always overwrite this setting using
  119 +
  120 +.. code-block:: bash
  121 +
  122 + export CALDB=/my/preferred/caldb/directory
  123 +
  124 +which you can add to your ``.bashrc`` file to set the directory
  125 +permanently.
  126 +ctools will use the ``CALDB`` environment variable to find out where the
  127 +calibration database is located. If needed, this mechanism can be
  128 +circumvented (see :ref:`sec_cta_rsp_abspath`).
  129 +
  130 +
  131 +Specifying the CTA Instrument Response Functions
  132 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  133 +
  134 +The specification of the CTA Instrument Response Functions depends on the
  135 +way how ctools are used. Common to all methods is that the IRFs are
  136 +defined by a response name and a calibration database name.
  137 +
  138 +There are different means to specify the CTA Instrument Response Functions
  139 +when using ctools, and the following section describe the
  140 +
  141 +
  142 +Specifying the response function as input parameters
  143 +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  144 +
  145 +ctools that require instrument response functions have two parameters
  146 +to specify the calibration database name and the response function name.
  147 +The following example shows a ``ctobssim`` run using the ``dummy``
  148 +calibration database and the ``cta_dummy_irf`` response function:
  149 +
  150 +.. code-block:: bash
  151 +
  152 + $ ctobssim
  153 + Model [$CTOOLS/share/models/crab.xml]
  154 + Calibration database [dummy]
  155 + Instrument response function [cta_dummy_irf]
  156 + RA of pointing (degrees) (0-360) [83.63]
  157 + Dec of pointing (degrees) (-90-90) [22.01]
  158 + Radius of FOV (degrees) (0-180) [5.0]
  159 + Start time (MET in s) (0) [0.0]
  160 + End time (MET in s) (0) [1800.0]
  161 + Lower energy limit (TeV) (0) [0.1]
  162 + Upper energy limit (TeV) (0) [100.0]
  163 + Output event data file or observation definition file [events.fits]
  164 +
  165 +Running the other tools is equivalent.
  166 +
  167 +
  168 +Specifying the response function in an observation definition file
  169 +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  170 +
  171 +In the above example, only a single global response function can be
  172 +used for all CTA observations. If you need to specify response functions
  173 +per observation you can add the information directly in the XML observation
  174 +definition file:
  175 +
  176 +.. code-block:: xml
  177 +
  178 + <observation_list title="observation library">
  179 + <observation name="Crab" id="00001" instrument="CTA">
  180 + <parameter name="EventList" file="events.fits"/>
  181 + <parameter name="Calibration" database="dummy" response="cta_dummy_irf"/>
  182 + </observation>
  183 + </observation_list>
  184 +
  185 +The ``Calibration`` parameter specifies the calibration database and
  186 +response name. You can then pass this file directly to, e.g., ``ctlike``:
  187 +
  188 +.. code-block:: bash
  189 +
  190 + $ ctlike
  191 + Event list, counts map or observation definition file [events.fits] obs_rsp.xml
  192 + Source model [$CTOOLS/share/models/crab.xml]
  193 + Source model output file [crab_results.xml]
  194 +
  195 +Note that ``ctlike`` does not ask for the calibration database and
  196 +response name as it found the relevant information in the XML file.
  197 +
  198 +.. _sec_cta_rsp_abspath:
  199 +
  200 +
  201 +Specifying individual instrument response files
  202 +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  203 +
  204 +If you need even more control over individual response files, you can
  205 +specify them individually in the XML observation file as follows:
  206 +
  207 +.. code-block:: xml
  208 +
  209 + <observation_list title="observation library">
  210 + <observation name="Crab" id="00001" instrument="CTA">
  211 + <parameter name="EventList" file="events.fits"/>
  212 + <parameter name="EffectiveArea" file="$CALDB/data/cta/dummy/bcf/cta_dummy_irf.dat"/>
  213 + <parameter name="PointSpreadFunction" file="$CALDB/data/cta/dummy/bcf/cta_dummy_irf.dat"/>
  214 + <parameter name="EnergyDispersion" file="$CALDB/data/cta/dummy/bcf/cta_dummy_irf.dat"/>
  215 + <parameter name="Background" file="$CALDB/data/cta/dummy/bcf/cta_dummy_irf.dat"/>
  216 + </observation>
  217 + </observation_list>
  218 +
  219 +
  220 +From within a Python script
  221 +^^^^^^^^^^^^^^^^^^^^^^^^^^^
  222 +
  223 +The following example illustrates how to set the calibration database
  224 +and response name from within Python:
  225 +
  226 +.. code-block:: python
  227 +
  228 + import gammalib
  229 + obs = gammalib.GCTAObservation()
  230 + caldb = gammalib.GCaldb("cta", "dummy")
  231 + irf = "cta_dummy_irf"
  232 + obs.response(irf, caldb)
  233 +
  234 +The calibration database is set by creating a ``GCaldb`` object. The
  235 +constructor takes as argument the mission (always ``cta``) and the
  236 +database name, in our case ``dummy``. The response function is then set
  237 +by passing the response name (here ``cta_dummy_irf``) and the calibration
  238 +database object to the ``response`` method.
  239 +
  240 +Alternatively, you can specify a specific response file by using
  241 +
  242 +.. code-block:: python
  243 +
  244 + import gammalib
  245 + obs = gammalib.GCTAObservation()
  246 + caldb = gammalib.GCaldb("$CALDB/data/cta/dummy/bcf")
  247 + irf = "cta_dummy_irf.dat"
  248 + obs.response(irf, caldb)
  249 +
  250 +Here, the directory in which the calibration file resides is specified
  251 +as argument to the ``GCaldb`` constructor, and the filename is passed
  252 +as response function name to the ``response`` method.
  253 +
  254 +
  255 +What response functions are available?
  256 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  257 +
  258 +To see what response functions are available on your system, type
  259 +
  260 +.. code-block:: bash
  261 +
  262 + $ cscaldb
  263 +
  264 +This creates the file ``cscaldb.log`` in your current working directory
  265 +that contains information about the available response functions. The
  266 +format of the file is:
  267 +
  268 +.. code-block:: bash
  269 +
  270 + 2014-02-20T17:15:21: +==============+
  271 + 2014-02-20T17:15:21: | Mission: cta |
  272 + 2014-02-20T17:15:21: +==============+
  273 + 2014-02-20T17:15:21: === Calibration: aar ===
  274 + 2014-02-20T17:15:21: DESY20140105_50h
  275 + 2014-02-20T17:15:21: DESY20140105_50h_0deg
  276 + 2014-02-20T17:15:21: DESY20140105_50h_180deg
  277 + 2014-02-20T17:15:21:
  278 + 2014-02-20T17:15:21: === Calibration: aar500 ===
  279 + 2014-02-20T17:15:21: DESY20140105_50h
  280 + 2014-02-20T17:15:21: DESY20140105_50h_0deg
  281 + 2014-02-20T17:15:21: DESY20140105_50h_180deg
  282 +
  283 +The name after ``Calibration`` is what we called earlier the calibration
  284 +database, and the names that are following each
  285 +``=== Calibration: xxx ===`` header are the response file names.
  286 +
  287 +To get the same information also onto the console, type
  288 +
  289 +.. code-block:: bash
  290 +
  291 + $ cscaldb debug=yes
  292 +
  293 +
  294 +Data format details
  295 +~~~~~~~~~~~~~~~~~~~
  296 +
  297 +.. _sec_cta_rsp_perftable:
  298 +
  299 +Performance tables
  300 +^^^^^^^^^^^^^^^^^^
  301 +
  302 +In the early days, the instrument performances derived from Monte-Carlo
  303 +simulations have been summarised in what we call here Performance
  304 +Tables, which are ASCII files that contain as function of energy the
  305 +on-axis performance parameters of CTA, such as effective area, point spread
  306 +function containment radius, energy resolution, background count rate and
  307 +differential sensitivity.
  308 +
  309 +Below an example of a CTA performance table::
  310 +
  311 + log(E) Area r68 r80 ERes. BG Rate Diff Sens
  312 + -1.7 261.6 0.3621 0.4908 0.5134 0.0189924 6.88237e-11
  313 + -1.5 5458.2 0.2712 0.3685 0.4129 0.1009715 1.72717e-11
  314 + -1.3 15590.0 0.1662 0.2103 0.2721 0.0575623 6.16963e-12
  315 + -1.1 26554.1 0.1253 0.1567 0.2611 0.0213008 2.89932e-12
  316 + -0.9 52100.5 0.1048 0.1305 0.1987 0.0088729 1.39764e-12
  317 + -0.7 66132.1 0.0827 0.1024 0.1698 0.0010976 6.03531e-13
  318 + -0.5 108656.8 0.0703 0.0867 0.1506 0.0004843 3.98147e-13
  319 + -0.3 129833.0 0.0585 0.0722 0.1338 0.0001575 3.23090e-13
  320 + -0.1 284604.3 0.0531 0.0656 0.1008 0.0001367 2.20178e-13
  321 + 0.1 263175.3 0.0410 0.0506 0.0831 0.0000210 1.87452e-13
  322 + 0.3 778048.6 0.0470 0.0591 0.0842 0.0000692 1.53976e-13
  323 + 0.5 929818.8 0.0391 0.0492 0.0650 0.0000146 1.18947e-13
  324 + 0.7 1078450.0 0.0335 0.0415 0.0541 0.0000116 1.51927e-13
  325 + 0.9 1448579.1 0.0317 0.0397 0.0516 0.0000047 1.42439e-13
  326 + 1.1 1899905.0 0.0290 0.0372 0.0501 0.0000081 1.96670e-13
  327 + 1.3 2476403.8 0.0285 0.0367 0.0538 0.0000059 2.20695e-13
  328 + 1.5 2832570.6 0.0284 0.0372 0.0636 0.0000073 3.22523e-13
  329 + 1.7 3534065.3 0.0290 0.0386 0.0731 0.0000135 4.84153e-13
  330 + 1.9 3250103.4 0.0238 0.0308 0.0729 0.0000044 6.26265e-13
  331 + 2.1 3916071.6 0.0260 0.0354 0.0908 0.0000023 7.69921e-13
  332 + ---------------------------------------------
  333 + 1) log(E) = log10(E/TeV) - bin centre
  334 + 2) Eff Area - in square metres after background cut (no theta cut)
  335 + 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees
  336 + 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees
  337 + 5) Fractional Energy Resolution (rms)
  338 + 6) BG Rate - inside point-source selection region - post call cuts - in Hz
  339 + 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE
  340 + - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including
  341 + systematics and statistics and at least 10 photons.
  342 +
  343 +
  344 +.. _sec_cta_rsp_xspec:
  345 +
  346 +ARF, RMF and PSF files
  347 +^^^^^^^^^^^^^^^^^^^^^^
  348 +
  349 +Instrument response information for the first CTA Data Challenge (1DC) has been
  350 +provided in a format that was heavily inspired by the
  351 +`RMF and ARF file formats
  352 +<http://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html>`_
  353 +that have been introduced by
  354 +`HEASARC <http://heasarc.gsfc.nasa.gov/>`_
  355 +for spectral analysis of X-ray data.
  356 +
  357 +The Redistribution Matrix File (RMF) describes how an incoming photon with
  358 +a given true energy is redistributed in measured energy. In other words, it
  359 +describes the energy dispersion of the instrument. The RMF is organised as
  360 +a two-dimensional matrix, with the first axis being in true energies while
  361 +the second axis represents the measured energies.
  362 +
  363 +The Ancillary Response File (ARF) describes the sensitivity of the
  364 +instrument to photons of a given true energy. The ARF gives the effective
  365 +area of the detector system after applying an event selection cut
  366 +(theta cut). For computing the ARF, some knowledge of the Point Spread
  367 +Function (PSF) is needed, so that the fraction of photons that falls into
  368 +the event selection region for a given source can be properly estimated.
  369 +Once this is done, no further PSF information is needed for the analysis.
  370 +
  371 +ctools however relies on the full modelling of the instrument, including
  372 +the Point Spread Function, and hence, PSF information has also been
  373 +provided for 1DC. Two data formats have been used for this: a first that
  374 +is based on a simple one-dimensional vector, providing the width of a
  375 +2-dimensional Gaussian as function of energy; and a second that is based
  376 +on a three-component 2D Gaussian as function of energy and offset angle
  377 +in the camera system. While the former was implement by a simply FITS
  378 +table with two columns, the latter was implemented by
  379 +:ref:`sec_cta_rsp_rsptable`.
  380 +
  381 +
  382 +.. _sec_cta_rsp_rsptable:
  383 +
  384 +Response tables
  385 +^^^^^^^^^^^^^^^
  386 +
  387 +The CTA response table class ``GCTAResponseTable`` provides a generic
  388 +handle for multi-dimensional response information. It is based on the
  389 +response format used for storing response information for the
  390 +*Fermi*/LAT telescope. In this format, all information is stored in
  391 +a single row of a FITS binary table. Each element of the row contains
  392 +a vector column, that describes the axes of the multi-dimensional response
  393 +cube and the response information. Note that this class may in the future
  394 +be promoted to the GammaLib core, as a similar class has been implemented
  395 +in the *Fermi*/LAT interface.
... ...
doc/source/user_manual/index.rst
  View file @89ab90b
... ... @@ -7,7 +7,7 @@ User Manual
7 7  
8 8 introduction
9 9 getting
10   - quickstart
  10 + getting_started/index
11 11 enhancing
12 12 reference
13 13 expert
... ...
scripts/Makefile.am
  View file @89ab90b
... ... @@ -4,12 +4,12 @@
4 4 binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py
5 5  
6 6 # Scripts to be distributed in pythondir
7   -pythonscripts = cspull.py cssens.py cstsdist.py \
  7 +pythonscripts = cspull.py cssens.py cstsdist.py cscaldb.py \
8 8 obsutils.py \
9 9 lsmodel.py
10 10  
11 11 # Scripts for which symbolic links will be created in bindir
12   -linkscripts = cspull cssens cstsdist \
  12 +linkscripts = cspull cssens cstsdist cscaldb \
13 13 lsmodel
14 14  
15 15 # Distribute scripts
... ...
scripts/cscaldb.py 0 → 100755
  View file @89ab90b
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script dumps all available calibrations into the console.
  4 +#
  5 +# Copyright (C) 2014 Juergen Knoedlseder
  6 +#
  7 +# This program is free software: you can redistribute it and/or modify
  8 +# it under the terms of the GNU General Public License as published by
  9 +# the Free Software Foundation, either version 3 of the License, or
  10 +# (at your option) any later version.
  11 +#
  12 +# This program is distributed in the hope that it will be useful,
  13 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 +# GNU General Public License for more details.
  16 +#
  17 +# You should have received a copy of the GNU General Public License
  18 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  19 +#
  20 +# ==========================================================================
  21 +import gammalib
  22 +import sys
  23 +import glob
  24 +import os
  25 +
  26 +
  27 +# ============ #
  28 +# cspull class #
  29 +# ============ #
  30 +class cscaldb(gammalib.GApplication):
  31 + """
  32 + This class implements the pull distribution generation script. It derives
  33 + from the GammaLib::GApplication class which provides support for parameter
  34 + files, command line arguments, and logging. In that way the Python
  35 + script behaves just as a regular ctool.
  36 + """
  37 + def __init__(self, *argv):
  38 + """
  39 + Constructor.
  40 + """
  41 + # Set name
  42 + self.name = "cscaldb"
  43 + self.version = "0.1.0"
  44 +
  45 + # Make sure that parfile exists
  46 + file = self.parfile()
  47 +
  48 + # Initialise application
  49 + if len(argv) == 0:
  50 + gammalib.GApplication.__init__(self, self.name, self.version)
  51 + elif len(argv) ==1:
  52 + gammalib.GApplication.__init__(self, self.name, self.version, *argv)
  53 + else:
  54 + raise TypeError("Invalid number of arguments given.")
  55 +
  56 + # Set logger properties
  57 + self.log_header()
  58 + self.log.date(True)
  59 +
  60 + # Return
  61 + return
  62 +
  63 + def __del__(self):
  64 + """
  65 + Destructor.
  66 + """
  67 + # Return
  68 + return
  69 +
  70 + def parfile(self):
  71 + """
  72 + Check if parfile exists. If parfile does not exist then create a
  73 + default parfile. This kluge avoids shipping the cscript with a parfile.
  74 + """
  75 + # Set parfile name
  76 + parfile = self.name+".par"
  77 +
  78 + try:
  79 + pars = gammalib.GApplicationPars(parfile)
  80 + except:
  81 + # Signal if parfile was not found
  82 + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
  83 +
  84 + # Create default parfile
  85 + pars = gammalib.GApplicationPars()
  86 + pars.append_standard()
  87 + pars.save(parfile)
  88 +
  89 + # Return
  90 + return
  91 +
  92 + def get_parameters(self):
  93 + """
  94 + Get parameters from parfile.
  95 + """
  96 + # Return
  97 + return
  98 +
  99 + def execute(self):
  100 + """
  101 + Execute the script.
  102 + """
  103 + # Run the script
  104 + self.run()
  105 +
  106 + # Return
  107 + return
  108 +
  109 + def run(self):
  110 + """
  111 + Run the script.
  112 + """
  113 + # Switch screen logging on in debug mode
  114 + if self.logDebug():
  115 + self.log.cout(True)
  116 +
  117 + # Get parameters
  118 + self.get_parameters()
  119 +
  120 + # Write input parameters into logger
  121 + if self.logTerse():
  122 + self.log_parameters()
  123 + self.log("\n")
  124 +
  125 + # Get the calibration database path
  126 + caldb = gammalib.GCaldb()
  127 +
  128 + # Gather missions
  129 + missions = []
  130 + paths = glob.glob(caldb.rootdir()+"/data/*")
  131 + paths.sort()
  132 + for path in paths:
  133 + head, tail = os.path.split(path)
  134 + missions.append(tail)
  135 +
  136 + # Loop over missions
  137 + for mission in missions:
  138 +
  139 + # Skip all non-CTA instruments
  140 + if mission != "cta":
  141 + continue
  142 +
  143 + # Write mission into logger
  144 + if self.logTerse():
  145 + self.log("\n")
  146 + self.log.header1("Mission: "+mission)
  147 +
  148 + # Gather instruments
  149 + instruments = []
  150 + paths = glob.glob(caldb.rootdir()+"/data/"+mission+"/*")
  151 + paths.sort()
  152 + for path in paths:
  153 + head, tail = os.path.split(path)
  154 + instruments.append(tail)
  155 +
  156 + # Loop over instruments
  157 + for instrument in instruments:
  158 +
  159 + # Write mission into logger
  160 + if self.logTerse():
  161 + self.log.header3("Calibration: "+instrument)
  162 +
  163 + # Open calibration index file
  164 + cifname = caldb.rootdir()+"/data/"+mission+"/"+instrument+"/caldb.indx"
  165 + fits = gammalib.GFits(cifname)
  166 + cif = fits["CIF"]
  167 + cals = cif["CAL_CBD"]
  168 +
  169 + # Extract response names
  170 + names = []
  171 + nrows = cals.length()
  172 + ncols = cals.number()
  173 + for row in range(nrows):
  174 + for col in range(ncols):
  175 + cal = cals.string(row, col)
  176 + istart = cal.find("NAME(")
  177 + if istart != -1:
  178 + istop = cal.find(")")
  179 + name = cal[5:istop]
  180 + if names.count(name) == 0:
  181 + names.append(name)
  182 + names.sort()
  183 +
  184 + # Print response name
  185 + if self.logTerse():
  186 + for name in names:
  187 + self.log(name+"\n")
  188 + self.log("\n")
  189 +
  190 + # Return
  191 + return
  192 +
  193 +
  194 +# ======================== #
  195 +# Main routine entry point #
  196 +# ======================== #
  197 +if __name__ == '__main__':
  198 + """
  199 + Dumps calibration database into log file.
  200 + """
  201 + # Create instance of application
  202 + app = cscaldb(sys.argv)
  203 +
  204 + # Open logfile
  205 + app.logFileOpen()
  206 +
  207 + # Execute application
  208 + app.execute()
  209 +
0 210 \ No newline at end of file
... ...
scripts/cssens.py
  View file @89ab90b
... ... @@ -57,7 +57,7 @@ class cssens(GApplication):
57 57 """
58 58 # Set name
59 59 self.name = "cssens"
60   - self.version = "0.3.0"
  60 + self.version = "0.4.0"
61 61  
62 62 # Make sure that parfile exists
63 63 file = self.parfile()
... ... @@ -106,11 +106,11 @@ class cssens(GApplication):
106 106 # Create default parfile
107 107 pars = GApplicationPars()
108 108 pars.append(GApplicationPar("outfile","f","h","sensitivity.dat","","","Output file name"))
109   - pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
  109 + pars.append(GApplicationPar("caldb","s","a","$CTOOLS/share/caldb/data/cta/dummy","","","Calibration database"))
110 110 pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
111 111 pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
112 112 pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
113   - pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
  113 + pars.append(GApplicationPar("bkg","s","a","$CTOOLS/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
114 114 pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
115 115 pars.append(GApplicationPar("rad","r","a","5.0","","","Radius of ROI (deg)"))
116 116 pars.append(GApplicationPar("enumbins","i","h","0","","","Number of energy bins (0=unbinned)"))
... ... @@ -333,7 +333,7 @@ class cssens(GApplication):
333 333 location.lb_deg(l, b)
334 334  
335 335 # Set source spectrum
336   - spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV")
  336 + spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV"))
337 337 if fitidx:
338 338 spectrum["Index"].free()
339 339 else:
... ...
scripts/cta_kb_root2caldb.py 0 → 100755
  View file @89ab90b
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script generates CALDB compliant performance tables from MC ROOT
  4 +# files.
  5 +#
  6 +# Note that the effective area given in the ASCII file is the full effective
  7 +# area (i.e. the area before applying any theta cut). This is required by
  8 +# GammaLib as the maximum likelihood analysis does not imply any theta cut.
  9 +#
  10 +# Copyright (C) 2014 Juergen Knoedlseder
  11 +#
  12 +# This program is free software: you can redistribute it and/or modify
  13 +# it under the terms of the GNU General Public License as published by
  14 +# the Free Software Foundation, either version 3 of the License, or
  15 +# (at your option) any later version.
  16 +#
  17 +# This program is distributed in the hope that it will be useful,
  18 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  19 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  20 +# GNU General Public License for more details.
  21 +#
  22 +# You should have received a copy of the GNU General Public License
  23 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  24 +#
  25 +# ==========================================================================
  26 +from ROOT import TFile, TH1F, TH2F
  27 +from datetime import datetime
  28 +import gammalib
  29 +import math
  30 +import sys
  31 +import glob
  32 +import os
  33 +
  34 +
  35 +# =========== #
  36 +# CalDB class #
  37 +# =========== #
  38 +class caldb():
  39 + """
  40 + """
  41 + def __init__(self, inst, path=""):
  42 + """
  43 + CALDB constructor. The constructor will create one CALDB entry for
  44 + an observation.
  45 +
  46 + Parameters:
  47 + inst - Instrument string.
  48 + Keywords:
  49 + path - CALDB Root path.
  50 + """
  51 + # Store root path
  52 + self.path = path
  53 +
  54 + # Initialise some members
  55 + self.cif = None
  56 + self.irf = None
  57 + self.ea = None
  58 + self.psf = None
  59 + self.edisp = None
  60 + self.bgd = None
  61 + self.hdu_cif = None
  62 +
  63 + # Store UTC date of job execution
  64 + self.utc = datetime.utcnow().isoformat()[:19]
  65 +
  66 + # Set CALDB information
  67 + self.cal_tel = "CTA"
  68 + self.cal_inst = inst.upper()
  69 + self.cal_det = "NONE"
  70 + self.cal_flt = "NONE"
  71 + self.cal_class = "BCF"
  72 + self.cal_type = "DATA"
  73 + self.cal_qual = 0 # 0=good, 1=bad, 2=dubious, ...
  74 + self.cal_date = "14/01/01"
  75 + self.val_date = "2014-01-01"
  76 + self.val_time = "00:00:00"
  77 + self.ref_time = 51544.0
  78 + self.cal_name = "NAME()"
  79 + self.cal_version = "VERSION(Prod1)"
  80 + self.cal_cut = "CLASS(BEST)"
  81 + self.cal_analysis = "ANALYSIS(MPIK)"
  82 + #
  83 + self.ea_name = "EFF_AREA"
  84 + self.ea_doc = "???"
  85 + self.ea_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
  86 + self.ea_desc = "CTA effective area"
  87 + #
  88 + self.psf_name = "RPSF"
  89 + self.psf_doc = "???"
  90 + self.psf_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
  91 + self.psf_desc = "CTA point spread function"
  92 + #
  93 + self.edisp_name = "EDISP"
  94 + self.edisp_doc = "???"
  95 + self.edisp_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
  96 + self.edisp_desc = "CTA energy dispersion"
  97 + #
  98 + self.bgd_name = "BGD"
  99 + self.bgd_doc = "???"
  100 + self.bgd_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
  101 + self.bgd_desc = "CTA background"
  102 +
  103 + # Create directory structure
  104 + self.make_dirs()
  105 +
  106 + # Return
  107 + return
  108 +
  109 + def __del__(self):
  110 + """
  111 + Destructor. Close all FITS files.
  112 + """
  113 + # Close calibration
  114 + self.close()
  115 +
  116 + # Return
  117 + return
  118 +
  119 + def make_dirs(self):
  120 + """
  121 + Generate CALDB directory structure for one observation identifier.
  122 + The structure is given by
  123 +
  124 + data/<tel>/<inst>/bcf/
  125 +
  126 + where <tel> is "cta" and <inst> is the instrument specified in the
  127 + CALDB constructor (the instrument may be used for different array
  128 + configurations).
  129 +
  130 + Parameters:
  131 + None
  132 + Keywords:
  133 + None
  134 + """
  135 + # Set base path
  136 + self.base_dir = "data"
  137 + self.base_dir += "/"+self.cal_tel.lower()
  138 + self.base_dir += "/"+self.cal_inst.lower()
  139 +
  140 + # Set directory names
  141 + self.bcf_dir = self.base_dir+"/bcf"
  142 +
  143 + # Optionally prefix path
  144 + if len(self.path) > 0:
  145 + self.base_path = self.path+"/"+self.base_dir
  146 + self.bcf_path = self.path+"/"+self.bcf_dir
  147 + else:
  148 + self.base_path = self.base_dir
  149 + self.bcf_path = self.bcf_dir
  150 +
  151 + # Create directory structure
  152 + if not os.path.isdir(self.bcf_path):
  153 + os.makedirs(self.bcf_path)
  154 +
  155 + # Return
  156 + return
  157 +
  158 + def add(self, rspname, split=False, clobber=True):
  159 + """
  160 + Add new calibration. The actual version will put
  161 + all calibrations in the same file, although each part of the response
  162 + function will have its own logical name. We can thus easily modify
  163 + the script to put each calibration information in a separate file.
  164 +
  165 + Parameters:
  166 + rspname - Response name
  167 + Keywords:
  168 + split - Split IRF over several files?
  169 + clobber - Overwrite existing files?
  170 + """
  171 + # Set calibrate file names
  172 + if split:
  173 + self.ea_file = "ea_"+rspname+".dat"
  174 + self.psf_file = "psf_"+rspname+".dat"
  175 + self.edisp_file = "edisp_"+rspname+".dat"
  176 + self.bgd_file = "bgd_"+rspname+".dat"
  177 + else:
  178 + self.ea_file = "irf_"+rspname+".dat"
  179 + self.psf_file = "irf_"+rspname+".dat"
  180 + self.edisp_file = "irf_"+rspname+".dat"
  181 + self.bgd_file = "irf_"+rspname+".dat"
  182 +
  183 + # Open calibration database index
  184 + if self.cif == None:
  185 + self.cif = gammalib.GFits(self.base_path+"/caldb.indx", True)
  186 +
  187 + # If file has no CIF extension than create it now
  188 + try:
  189 + self.hdu_cif = self.cif.table("CIF")
  190 + except:
  191 + self.create_cif()
  192 + self.hdu_cif = self.cif.table("CIF")
  193 +
  194 + # Set response name
  195 + self.cal_name = "NAME("+rspname+")"
  196 +
  197 + # Return
  198 + return
  199 +
  200 + def close(self):
  201 + """
  202 + Close calibration FITS files.
  203 +
  204 + Parameters:
  205 + None
  206 + Keywords:
  207 + None
  208 + """
  209 + # Add information to CIF. We do this now as before we did not
  210 + # necessarily have all the information at hand (in particular
  211 + # about the boundaries)
  212 + #self.add_cif_info()
  213 +
  214 + # Close CIF
  215 + if self.cif != None:
  216 + self.cif.save(True)
  217 + self.cif.close()
  218 + self.cif = None
  219 +
  220 + # Free
  221 + self.cif = None
  222 +
  223 + # Return
  224 + return
  225 +
  226 + def create_cif(self):
  227 + """
  228 + Create Calibration Index File (CIF) extension in CIF FITS file.
  229 +
  230 + Parameters:
  231 + None
  232 + Keywords:
  233 + None
  234 + """
  235 + # Create binary table
  236 + table = gammalib.GFitsBinTable()
  237 +
  238 + # Set boundary
  239 +
  240 + # Attach columns. Reference: CAL/GEN/92-008
  241 + table.append(gammalib.GFitsTableStringCol("TELESCOP", 0, 10))
  242 + table.append(gammalib.GFitsTableStringCol("INSTRUME", 0, 10))
  243 + table.append(gammalib.GFitsTableStringCol("DETNAM", 0, 20))
  244 + table.append(gammalib.GFitsTableStringCol("FILTER", 0, 10))
  245 + table.append(gammalib.GFitsTableStringCol("CAL_DEV", 0, 20))
  246 + table.append(gammalib.GFitsTableStringCol("CAL_DIR", 0, 70))
  247 + table.append(gammalib.GFitsTableStringCol("CAL_FILE", 0, 40))
  248 + table.append(gammalib.GFitsTableStringCol("CAL_CLAS", 0, 3))
  249 + table.append(gammalib.GFitsTableStringCol("CAL_DTYP", 0, 4))
  250 + table.append(gammalib.GFitsTableStringCol("CAL_CNAM", 0, 20))
  251 + table.append(gammalib.GFitsTableStringCol("CAL_CBD", 0, 70, 9))
  252 + table.append(gammalib.GFitsTableShortCol("CAL_XNO", 0))
  253 + table.append(gammalib.GFitsTableStringCol("CAL_VSD", 0, 10))
  254 + table.append(gammalib.GFitsTableStringCol("CAL_VST", 0, 8))
  255 + table.append(gammalib.GFitsTableDoubleCol("REF_TIME", 0))
  256 + table.append(gammalib.GFitsTableShortCol("CAL_QUAL", 0))
  257 + table.append(gammalib.GFitsTableStringCol("CAL_DATE", 0, 8))
  258 + table.append(gammalib.GFitsTableStringCol("CAL_DESC", 0, 70))
  259 +
  260 + # Set keywords. Reference: CAL/GEN/92-008
  261 + table.extname("CIF")
  262 + table.card("CIFVERSN", "1992a", "Version of CIF format")
  263 +
  264 + # Attach table to FITS file. Note that at this moment the
  265 + # FITS table gets connected to the FITS file. Yet since
  266 + # nothing was yet written to the FITS file, we cannot read
  267 + # anything from it.
  268 + self.cif.append(table)
  269 +
  270 + # Return
  271 + return
  272 +
  273 + def add_cif_info(self):
  274 + """
  275 + Add information to CIF extension.
  276 +
  277 + Parameters:
  278 + None
  279 + Keywords:
  280 + None
  281 + """
  282 + # Set boundaries
  283 + self.ea_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
  284 + self.psf_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
  285 + self.edisp_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
  286 + self.bgd_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
  287 +
  288 + # Append 4 rows to CIF extension
  289 + self.hdu_cif.append_rows(4)
  290 +
  291 + # Add generic information for these 4 rows
  292 + for i in range(4):
  293 +
  294 + # Set row number
  295 + row = i+self.hdu_cif.nrows()-4
  296 +
  297 + # Set element
  298 + self.hdu_cif["TELESCOP"][row] = self.cal_tel
  299 + self.hdu_cif["INSTRUME"][row] = self.cal_inst
  300 + self.hdu_cif["DETNAM"][row] = self.cal_det
  301 + self.hdu_cif["FILTER"][row] = self.cal_flt
  302 + self.hdu_cif["CAL_DEV"][row] = "ONLINE"
  303 + self.hdu_cif["CAL_CLAS"][row] = self.cal_class
  304 + self.hdu_cif["CAL_DTYP"][row] = self.cal_type
  305 + self.hdu_cif["CAL_VSD"][row] = self.val_date
  306 + self.hdu_cif["CAL_VST"][row] = self.val_time
  307 + self.hdu_cif["REF_TIME"][row] = self.ref_time
  308 + self.hdu_cif["CAL_QUAL"][row] = self.cal_qual
  309 + self.hdu_cif["CAL_DATE"][row] = self.cal_date
  310 +
  311 + # Add effective area information
  312 + row = self.hdu_cif.nrows()-4
  313 + self.hdu_cif["CAL_DIR"][row] = self.bcf_dir
  314 + self.hdu_cif["CAL_FILE"][row] = self.ea_file
  315 + self.hdu_cif["CAL_CNAM"][row] = self.ea_name
  316 + self.hdu_cif["CAL_DESC"][row] = self.ea_desc
  317 + self.hdu_cif["CAL_XNO"][row] = 1
  318 + for n in range(9):
  319 + if n >= len(self.ea_bounds):
  320 + self.hdu_cif["CAL_CBD"][row, n] = "NONE"
  321 + else:
  322 + self.hdu_cif["CAL_CBD"][row, n] = self.ea_bounds[n]
  323 +
  324 + # Add point spread function information
  325 + row = self.hdu_cif.nrows()-3
  326 + self.hdu_cif["CAL_DIR"][row] = self.bcf_dir
  327 + self.hdu_cif["CAL_FILE"][row] = self.psf_file
  328 + self.hdu_cif["CAL_CNAM"][row] = self.psf_name
  329 + self.hdu_cif["CAL_DESC"][row] = self.psf_desc
  330 + self.hdu_cif["CAL_XNO"][row] = 1
  331 + for n in range(9):
  332 + if n >= len(self.psf_bounds):
  333 + self.hdu_cif["CAL_CBD"][row, n] = "NONE"
  334 + else:
  335 + self.hdu_cif["CAL_CBD"][row, n] = self.psf_bounds[n]
  336 +
  337 + # Add energy dispersion information
  338 + row = self.hdu_cif.nrows()-2
  339 + self.hdu_cif["CAL_DIR"][row] = self.bcf_dir
  340 + self.hdu_cif["CAL_FILE"][row] = self.edisp_file
  341 + self.hdu_cif["CAL_CNAM"][row] = self.edisp_name
  342 + self.hdu_cif["CAL_DESC"][row] = self.edisp_desc
  343 + self.hdu_cif["CAL_XNO"][row] = 1
  344 + for n in range(9):
  345 + if n >= len(self.edisp_bounds):
  346 + self.hdu_cif["CAL_CBD"][row, n] = "NONE"
  347 + else:
  348 + self.hdu_cif["CAL_CBD"][row, n] = self.edisp_bounds[n]
  349 +
  350 + # Add background information
  351 + row = self.hdu_cif.nrows()-1
  352 + self.hdu_cif["CAL_DIR"][row] = self.bcf_dir
  353 + self.hdu_cif["CAL_FILE"][row] = self.bgd_file
  354 + self.hdu_cif["CAL_CNAM"][row] = self.bgd_name
  355 + self.hdu_cif["CAL_DESC"][row] = self.bgd_desc
  356 + self.hdu_cif["CAL_XNO"][row] = 1
  357 + for n in range(9):
  358 + if n >= len(self.bgd_bounds):
  359 + self.hdu_cif["CAL_CBD"][row, n] = "NONE"
  360 + else:
  361 + self.hdu_cif["CAL_CBD"][row, n] = self.bgd_bounds[n]
  362 +
  363 + # Return
  364 + return
  365 +
  366 + def root2caldb(self, filename):
  367 + """
  368 + Generate files.
  369 +
  370 + Parameters:
  371 + filename - Performance file name
  372 + """
  373 + # Open ROOT performance file
  374 + rootfile = TFile(filename)
  375 +
  376 + # Set filenames
  377 + ea_filename = self.bcf_path+"/"+self.ea_file
  378 + psf_filename = self.bcf_path+"/"+self.psf_file
  379 + edisp_filename = self.bcf_path+"/"+self.edisp_file
  380 + bgd_filename = self.bcf_path+"/"+self.bgd_file
  381 +
  382 + # Open ASCII performance tables
  383 + if self.ea_file == self.psf_file and \
  384 + self.ea_file == self.edisp_file and \
  385 + self.ea_file == self.bgd_file :
  386 + split = False
  387 + fprf = open(ea_filename, "w")
  388 + files = [fprf]
  389 + else:
  390 + split = True
  391 + fea = open(ea_filename, "w")
  392 + fpsf = open(psf_filename, "w")
  393 + fedisp = open(edisp_filename, "w")
  394 + fbgd = open(bgd_filename, "w")
  395 + files = [fea, fpsf, fedisp, fbgd]
  396 +
  397 + # Write header
  398 + for file in files:
  399 + file.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n")
  400 +
  401 + # Get histograms. Konrad's root files do not have an EffectiveArea80
  402 + # histogram, but his effective area is for 80% containment radius
  403 + # anyways, so we simply read this histogram into aeff80.
  404 + sens = TH1F()
  405 + bgrate = TH1F()
  406 + bgratesqdeg = TH1F()
  407 + aeff = TH1F()
  408 + aeff80 = TH1F()
  409 + angres = TH1F()
  410 + angres80 = TH1F()
  411 + eres = TH1F()
  412 + rootfile.GetObject("DiffSens", sens)
  413 + rootfile.GetObject("BGRate", bgrate)
  414 + rootfile.GetObject("BGRatePerSqDeg", bgratesqdeg)
  415 + rootfile.GetObject("EffectiveArea", aeff)
  416 + try:
  417 + rootfile.GetObject("EffectiveArea80", aeff80)
  418 + except:
  419 + rootfile.GetObject("EffectiveArea", aeff80)
  420 + rootfile.GetObject("AngRes", angres)
  421 + rootfile.GetObject("AngRes80", angres80)
  422 + rootfile.GetObject("ERes", eres)
  423 +
  424 + # Setup vectors
  425 + energies = bgrate.GetXaxis()
  426 + nbins_eng = energies.GetNbins()
  427 + for i in range(nbins_eng):
  428 +
  429 + # Get logE and background rate per squared degree
  430 + logE = energies.GetBinCenter(i+1)
  431 + bgd = bgrate.GetBinContent(i+1,1)
  432 + bgdsq = bgratesqdeg.GetBinContent(i+1,1)
  433 + area = aeff.GetBinContent(i+1,1)
  434 + area80 = aeff80.GetBinContent(i+1,1)
  435 + r68 = angres.GetBinContent(i+1,1)
  436 + r80 = angres80.GetBinContent(i+1,1)
  437 + dE = eres.GetBinContent(i+1,1)
  438 + diffs = sens.GetBinContent(i+1,1)
  439 +
  440 + # Skip any NaNs
  441 + if math.isnan(area80):
  442 + continue
  443 +
  444 + # Compute energy (in MeV)
  445 + energy = math.pow(10.0, logE)*1.0e6
  446 + emin = math.pow(10.0, logE-0.1)*1.0e6
  447 + emax = math.pow(10.0, logE+0.1)*1.0e6
  448 + ewidth = emax-emin
  449 +
  450 + # Convert into background rate per steradian and MeV
  451 + bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth
  452 +
  453 + # Compute control background rate (this only works for KB files
  454 + # as they are for 80% containment radius). But we don't use this
  455 + # anyways, it's just for control on display ...
  456 + omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80)))
  457 + if omega > 0.0:
  458 + bkg_rate2 = bgd / omega / ewidth
  459 + else:
  460 + bkg_rate2 = 0.0
  461 +
  462 + # Compute full effective area by dividing area80/0.80
  463 + area = area80/0.80
  464 +
  465 + # Write results in file
  466 + line = "%.1f %.1f %.4f %.4f %.4f %.5e %.5e" % \
  467 + (logE, area, r68, r80, dE, bgd, diffs)
  468 + for file in files:
  469 + file.write(line+"\n")
  470 + #fbkg.write(str(energy)+" "+str(bkg_rate)+"\n")
  471 +
  472 + # Write trailer
  473 + for file in files:
  474 + file.write("---------------------------------------------\n")
  475 + file.write("Notes\n")
  476 + file.write(" 1) log(E) = log10(E/TeV) - bin centre\n")
  477 + file.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n")
  478 + file.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n")
  479 + file.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n")
  480 + file.write(" 5) Fractional Energy Resolution (rms)\n")
  481 + file.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n")
  482 + file.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n")
  483 + file.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n")
  484 + file.write(" systematics and statistics and at least 10 photons.\n")
  485 +
  486 + # Close files
  487 + for file in files:
  488 + file.close()
  489 +
  490 + # Close ROOT performance file
  491 + rootfile.Close()
  492 +
  493 + # Add information to CIF. We do this now as before we did not
  494 + # necessarily have all the information at hand (in particular
  495 + # about the boundaries)
  496 + self.add_cif_info()
  497 +
  498 + # Return
  499 + return
  500 +
  501 +
  502 +#==========================#
  503 +# Main routine entry point #
  504 +#==========================#
  505 +if __name__ == '__main__':
  506 + """
  507 + """
  508 + # Get ROOT filename from command line argument
  509 + if (len(sys.argv) > 1):
  510 + filename = sys.argv[1]
  511 + else:
  512 + db = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_20deg_root"
  513 + #db = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_50deg_root"
  514 +
  515 + # Get list of all performance files
  516 + files = glob.glob(db+"/*.root")
  517 + files.sort()
  518 +
  519 + # Allocate caldb
  520 + irf = caldb("kb")
  521 +
  522 + # Loop over all files
  523 + for file in files:
  524 +
  525 + # Get filename
  526 + head, tail = os.path.split(file)
  527 +
  528 + # Build name
  529 + name = tail.strip(".root")
  530 + if name.find("kb_") != -1:
  531 + name = name.replace("_20deg", "")
  532 + elif name.find("IFAE_") != -1:
  533 + name = name.replace("Subarray", "")
  534 + name = name.replace("IFAE_", "")
  535 + name = "IFAE_"+name
  536 + i = name.find("hours")
  537 + if i != -1:
  538 + name = name[0:i+1]
  539 + print(tail)
  540 +
  541 + # Open calibration file
  542 + irf.add(name)
  543 +
  544 + # Translate ROOT to CALDB information
  545 + irf.root2caldb(file)
... ...
scripts/cta_root2caldb.py
  View file @89ab90b
... ... @@ -1319,9 +1319,9 @@ def set_prod1_ifae():
1319 1319 #=========================#
1320 1320 # Set Prod2 DESY database #
1321 1321 #=========================#
1322   -def set_prod2_desy():
  1322 +def set_prod2_desy_ignore():
1323 1323 """
1324   - Set Prod2 DESY database.
  1324 + Set Prod2 DESY database (ignore these files).
1325 1325 """
1326 1326 # Set database attributes
1327 1327 path_aar = "/project-data/cta/performance/prod2/Performance_DESY_20140131_Aar"
... ... @@ -1366,6 +1366,56 @@ def set_prod2_desy():
1366 1366 return db
1367 1367  
1368 1368  
  1369 +#=========================#
  1370 +# Set Prod2 DESY database #
  1371 +#=========================#
  1372 +def set_prod2_desy():
  1373 + """
  1374 + Set Prod2 DESY database.
  1375 + """
  1376 + # Set database attributes
  1377 + path_aar = "/project-data/cta/performance/prod2/Performance_DESY_20140131_Aar"
  1378 + path_ten = "/project-data/cta/performance/prod2/Performance_DESY_20140131_Tenerife_50h"
  1379 + rebin = False
  1380 + psftype = "Gauss"
  1381 + scale = 1.0
  1382 +
  1383 + # Set database content
  1384 + db = [{'inst': "aar", 'id': "DESY20140105_50h",
  1385 + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1386 + 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-NS.S.2a.180000s.root"},
  1387 + {'inst': "aar", 'id': "DESY20140105_50h_0deg",
  1388 + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1389 + 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
  1390 + {'inst': "aar", 'id': "DESY20140105_50h_180deg",
  1391 + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1392 + 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
  1393 +
  1394 + {'inst': "aar500", 'id': "DESY20140105_50h",
  1395 + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1396 + 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-500m-NS.S.2a.180000s.root"},
  1397 + {'inst': "aar500", 'id': "DESY20140105_50h_0deg",
  1398 + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1399 + 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
  1400 + {'inst': "aar500", 'id': "DESY20140105_50h_180deg",
  1401 + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1402 + 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
  1403 +
  1404 + {'inst': "tenerife", 'id': "DESY20140105_50h",
  1405 + 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1406 + 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Tenerife-NS.N.2NN.180000s.root"},
  1407 + {'inst': "tenerife", 'id': "DESY20140105_50h_0deg",
  1408 + 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1409 + 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"},
  1410 + {'inst': "tenerife", 'id': "DESY20140105_50h_180deg",
  1411 + 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1412 + 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"}
  1413 + ]
  1414 +
  1415 + # Return database
  1416 + return db
  1417 +
  1418 +
1369 1419 #==========================#
1370 1420 # Main routine entry point #
1371 1421 #==========================#
... ... @@ -1381,7 +1431,7 @@ if __name__ == &#39;__main__&#39;:
1381 1431 - EffectiveArea80_offaxis (2D) - Effective area for 80% source containment
1382 1432 - AngRes_offaxis (2D) - Angular resolution 68% containment
1383 1433 - AngRes80_offaxis (2D) - Angular resolution 80% containment
1384   - - BGRatePerSqDeg_offaxis (2D) - Background rate per square degree
  1434 + - BGRatePerSqDeg_offaxis (2D) - Background rate per square degree [counts/s/deg2]
1385 1435 - BGRate_offaxis (2D) - Background rate
1386 1436 - EffectiveAreaEtrue_offaxis (2D) - Effective area in true energy (finer bins)
1387 1437 - MigMatrix_offaxis (3D)
... ...
scripts/cta_root2irf.py
  View file @89ab90b
... ... @@ -38,58 +38,58 @@ import os
38 38 # Compute Eq. (17) of Li & Ma (1983) #
39 39 # ================================== #
40 40 def sigma_lima(Non, Noff, alpha=0.2):
41   - """
42   - Compute Eq. (17) of Li & Ma (1983).
43   -
44   - Parameters:
45   - Non - Number of on counts
46   - Noff - Number of off counts
47   - Keywords:
48   - alpha - Ratio of on-to-off exposure
49   - """
50   - # Compute sensitivity
51   - alpha1 = alpha + 1.0
52   - sum = Non + Noff
53   - arg1 = Non / sum
54   - arg2 = Noff / sum
55   - term1 = Non * math.log((alpha1/alpha)*arg1)
56   - term2 = Noff * math.log(alpha1*arg2)
57   - sigma = math.sqrt(2.0 * (term1 + term2))
58   -
59   - # Return sensitivity
60   - return sigma
  41 + """
  42 + Compute Eq. (17) of Li & Ma (1983).
  43 +
  44 + Parameters:
  45 + Non - Number of on counts
  46 + Noff - Number of off counts
  47 + Keywords:
  48 + alpha - Ratio of on-to-off exposure
  49 + """
  50 + # Compute sensitivity
  51 + alpha1 = alpha + 1.0
  52 + sum = Non + Noff
  53 + arg1 = Non / sum
  54 + arg2 = Noff / sum
  55 + term1 = Non * math.log((alpha1/alpha)*arg1)
  56 + term2 = Noff * math.log(alpha1*arg2)
  57 + sigma = math.sqrt(2.0 * (term1 + term2))
  58 +
  59 + # Return sensitivity
  60 + return sigma
61 61  
62 62  
63 63 # ======================================== #
64 64 # Solve Eq. (17) of Li & Ma (1983) for Non #
65 65 # ======================================== #
66 66 def Non_lima(sigma, Noff, alpha=0.2):
67   - """
68   - Solve Eq. (17) of Li & Ma (1983) for Non
69   -
70   - Parameters:
71   - sigma - Required significance
72   - Noff - Number of off counts
73   - Keywords:
74   - alpha - Ratio of on-to-off exposure
75   - """
76   - # Set initial guess for Non
77   - S = 25.0
78   - Non = S + alpha*Noff
79   -
80   - # Compute initial sigma estimate
81   - s = sigma_lima(Non, Noff, alpha=alpha)
82   - delta = s - sigma
83   -
84   - # Iterate
85   - while abs(delta) > 1.0e-6:
86   - S *= sigma/s
87   - Non = S + alpha*Noff
88   - s = sigma_lima(Non, Noff, alpha=alpha)
89   - delta = s - sigma
90   -
91   - # Return
92   - return Non
  67 + """
  68 + Solve Eq. (17) of Li & Ma (1983) for Non
  69 +
  70 + Parameters:
  71 + sigma - Required significance
  72 + Noff - Number of off counts
  73 + Keywords:
  74 + alpha - Ratio of on-to-off exposure
  75 + """
  76 + # Set initial guess for Non
  77 + S = 25.0
  78 + Non = S + alpha*Noff
  79 +
  80 + # Compute initial sigma estimate
  81 + s = sigma_lima(Non, Noff, alpha=alpha)
  82 + delta = s - sigma
  83 +
  84 + # Iterate
  85 + while abs(delta) > 1.0e-6:
  86 + S *= sigma/s
  87 + Non = S + alpha*Noff
  88 + s = sigma_lima(Non, Noff, alpha=alpha)
  89 + delta = s - sigma
  90 +
  91 + # Return
  92 + return Non
93 93  
94 94  
95 95 # =============================================== #
... ... @@ -97,195 +97,195 @@ def Non_lima(sigma, Noff, alpha=0.2):
97 97 # =============================================== #
98 98 def get_aeff(bgrate, emin, emax, \
99 99 duration=180000.0, alpha=0.2, sigma=5.0, norm=2.79e-7, index=-2.57):
100   - """
101   - Computes effective area assuming a Crab spectrum. This does not
102   - include any hadron cut efficiencies, and thus is not very useful in a first
103   - place ...
104   -
105   - Parameters:
106   - bgrate - Background rate (TeV)
107   - emin - Minimum energy (TeV)
108   - emax - Maximum energy (TeV)
109   - Keywords:
110   - duration - Observing time in seconds (default: 50h)
111   - alpha - Non/Noff ratio (default: 0.2)
112   - sigma - Significance (default: 5.0)
113   - norm - Crab spectrum normalization at 1 TeV (default: 2.79e-7)
114   - index - Crab spectral index (default: -2.57)
115   - """
116   - # Compute number of source events
117   - B = bgrate*duration
118   - Noff = B/alpha
119   - Non = Non_lima(sigma, Noff)
120   - S = (Non - B)
121   -
122   - # Compute Crab flux in specified energy band
123   - inx = index+1.0
124   - cu = norm/inx * (math.pow(emax, inx) - math.pow(emin, inx))
125   -
126   - # Compute effective area
127   - area = S / duration / cu
128   -
129   - # Return effective area
130   - return area
  100 + """
  101 + Computes effective area assuming a Crab spectrum. This does not
  102 + include any hadron cut efficiencies, and thus is not very useful in a first
  103 + place ...
  104 +
  105 + Parameters:
  106 + bgrate - Background rate (TeV)
  107 + emin - Minimum energy (TeV)
  108 + emax - Maximum energy (TeV)
  109 + Keywords:
  110 + duration - Observing time in seconds (default: 50h)
  111 + alpha - Non/Noff ratio (default: 0.2)
  112 + sigma - Significance (default: 5.0)
  113 + norm - Crab spectrum normalization at 1 TeV (default: 2.79e-7)
  114 + index - Crab spectral index (default: -2.57)
  115 + """
  116 + # Compute number of source events
  117 + B = bgrate*duration
  118 + Noff = B/alpha
  119 + Non = Non_lima(sigma, Noff)
  120 + S = (Non - B)
  121 +
  122 + # Compute Crab flux in specified energy band
  123 + inx = index+1.0
  124 + cu = norm/inx * (math.pow(emax, inx) - math.pow(emin, inx))
  125 +
  126 + # Compute effective area
  127 + area = S / duration / cu
  128 +
  129 + # Return effective area
  130 + return area
131 131  
132 132  
133 133 # ============== #
134 134 # Generate files #
135 135 # ============== #
136 136 def make_files(filename, prfname, bkgname):
137   - """
138   - Generate files.
139   -
140   - Parameters:
141   - filename - Performance file name
142   - prfname - Performance file name
143   - bkgname - File function filename
144   - """
145   - # Open performance file and file function
146   - file = TFile(filename)
147   - fprf = open(prfname, "w")
148   - fbkg = open(bkgname, "w")
149   -
150   - # Write header
151   - fprf.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n")
152   - #sys.stdout.write("log(E) Area Area80 r68 r80 ERes. BG Rate Diff Sens BGRateSqDeg BGInt (BGIntControl)\n")
  137 + """
  138 + Generate files.
  139 +
  140 + Parameters:
  141 + filename - Performance file name
  142 + prfname - Performance file name
  143 + bkgname - File function filename
  144 + """
  145 + # Open performance file and file function
  146 + file = TFile(filename)
  147 + fprf = open(prfname, "w")
  148 + fbkg = open(bkgname, "w")
  149 +
  150 + # Write header
  151 + fprf.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n")
  152 + #sys.stdout.write("log(E) Area Area80 r68 r80 ERes. BG Rate Diff Sens BGRateSqDeg BGInt (BGIntControl)\n")
153 153  
154   - # Get histograms. Konrad's root files do not have an EffectiveArea80
155   - # histogram, but his effective area is for 80% containment radius
156   - # anyways, so we simply read this histogram into aeff80.
157   - sens = TH1F()
158   - bgrate = TH1F()
159   - bgratesqdeg = TH1F()
160   - aeff = TH1F()
161   - aeff80 = TH1F()
162   - angres = TH1F()
163   - angres80 = TH1F()
164   - eres = TH1F()
165   - file.GetObject("DiffSens", sens)
166   - file.GetObject("BGRate", bgrate)
167   - file.GetObject("BGRatePerSqDeg", bgratesqdeg)
168   - file.GetObject("EffectiveArea", aeff)
169   - try:
170   - file.GetObject("EffectiveArea80", aeff80)
171   - except:
172   - file.GetObject("EffectiveArea", aeff80)
173   - file.GetObject("AngRes", angres)
174   - file.GetObject("AngRes80", angres80)
175   - file.GetObject("ERes", eres)
176   -
177   - # Setup vectors
178   - energies = bgrate.GetXaxis()
179   - nbins_eng = energies.GetNbins()
180   - for i in range(nbins_eng):
181   -
182   - # Get logE and background rate per squared degree
183   - logE = energies.GetBinCenter(i+1)
184   - bgd = bgrate.GetBinContent(i+1,1)
185   - bgdsq = bgratesqdeg.GetBinContent(i+1,1)
186   - area = aeff.GetBinContent(i+1,1)
187   - area80 = aeff80.GetBinContent(i+1,1)
188   - r68 = angres.GetBinContent(i+1,1)
189   - r80 = angres80.GetBinContent(i+1,1)
190   - dE = eres.GetBinContent(i+1,1)
191   - diffs = sens.GetBinContent(i+1,1)
192   -
193   - # Skip any NaNs
194   - if math.isnan(area80):
195   - continue
  154 + # Get histograms. Konrad's root files do not have an EffectiveArea80
  155 + # histogram, but his effective area is for 80% containment radius
  156 + # anyways, so we simply read this histogram into aeff80.
  157 + sens = TH1F()
  158 + bgrate = TH1F()
  159 + bgratesqdeg = TH1F()
  160 + aeff = TH1F()
  161 + aeff80 = TH1F()
  162 + angres = TH1F()
  163 + angres80 = TH1F()
  164 + eres = TH1F()
  165 + file.GetObject("DiffSens", sens)
  166 + file.GetObject("BGRate", bgrate)
  167 + file.GetObject("BGRatePerSqDeg", bgratesqdeg)
  168 + file.GetObject("EffectiveArea", aeff)
  169 + try:
  170 + file.GetObject("EffectiveArea80", aeff80)
  171 + except:
  172 + file.GetObject("EffectiveArea", aeff80)
  173 + file.GetObject("AngRes", angres)
  174 + file.GetObject("AngRes80", angres80)
  175 + file.GetObject("ERes", eres)
  176 +
  177 + # Setup vectors
  178 + energies = bgrate.GetXaxis()
  179 + nbins_eng = energies.GetNbins()
  180 + for i in range(nbins_eng):
  181 +
  182 + # Get logE and background rate per squared degree
  183 + logE = energies.GetBinCenter(i+1)
  184 + bgd = bgrate.GetBinContent(i+1,1)
  185 + bgdsq = bgratesqdeg.GetBinContent(i+1,1)
  186 + area = aeff.GetBinContent(i+1,1)
  187 + area80 = aeff80.GetBinContent(i+1,1)
  188 + r68 = angres.GetBinContent(i+1,1)
  189 + r80 = angres80.GetBinContent(i+1,1)
  190 + dE = eres.GetBinContent(i+1,1)
  191 + diffs = sens.GetBinContent(i+1,1)
  192 +
  193 + # Skip any NaNs
  194 + if math.isnan(area80):
  195 + continue
196 196  
197   - # Compute energy (in MeV)
198   - energy = math.pow(10.0, logE)*1.0e6
199   - emin = math.pow(10.0, logE-0.1)*1.0e6
200   - emax = math.pow(10.0, logE+0.1)*1.0e6
201   - ewidth = emax-emin
202   -
203   - # Convert into background rate per steradian and MeV
204   - bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth
205   -
206   - # Compute control background rate (this only works for KB files
207   - # as they are for 80% containment radius). But we don't use this
208   - # anyways, it's just for control on display ...
209   - omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80)))
210   - if omega > 0.0:
211   - bkg_rate2 = bgd / omega / ewidth
212   - else:
213   - bkg_rate2 = 0.0
214   -
215   - # Compute full effective area by dividing area80/0.80
216   - area = area80/0.80
  197 + # Compute energy (in MeV)
  198 + energy = math.pow(10.0, logE)*1.0e6
  199 + emin = math.pow(10.0, logE-0.1)*1.0e6
  200 + emax = math.pow(10.0, logE+0.1)*1.0e6
  201 + ewidth = emax-emin
  202 +
  203 + # Convert into background rate per steradian and MeV
  204 + bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth
  205 +
  206 + # Compute control background rate (this only works for KB files
  207 + # as they are for 80% containment radius). But we don't use this
  208 + # anyways, it's just for control on display ...
  209 + omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80)))
  210 + if omega > 0.0:
  211 + bkg_rate2 = bgd / omega / ewidth
  212 + else:
  213 + bkg_rate2 = 0.0
  214 +
  215 + # Compute full effective area by dividing area80/0.80
  216 + area = area80/0.80
217 217  
218   - # Write results in file
219   - line = "%.1f %.1f %.4f %.4f %.4f %.5e %.5e" % \
220   - (logE, area, r68, r80, dE, bgd, diffs)
221   - fprf.write(line+"\n")
222   - fbkg.write(str(energy)+" "+str(bkg_rate)+"\n")
  218 + # Write results in file
  219 + line = "%.1f %.1f %.4f %.4f %.4f %.5e %.5e" % \
  220 + (logE, area, r68, r80, dE, bgd, diffs)
  221 + fprf.write(line+"\n")
  222 + fbkg.write(str(energy)+" "+str(bkg_rate)+"\n")
223 223  
224   - # Show performance file
225   - #print "%.1f %.1f %.1f %.4f %.4f %.4f %.7f %.5e %.5e %.5e (%.5e)" % \
226   - # (logE, area, area80, r68, r80, dE, bgd, diffs, bgdsq, bkg_rate, bkg_rate2)
  224 + # Show performance file
  225 + #print "%.1f %.1f %.1f %.4f %.4f %.4f %.7f %.5e %.5e %.5e (%.5e)" % \
  226 + # (logE, area, area80, r68, r80, dE, bgd, diffs, bgdsq, bkg_rate, bkg_rate2)
227 227  
228   - # Write trailer
229   - fprf.write("---------------------------------------------\n")
230   - fprf.write("Notes\n")
231   - fprf.write(" 1) log(E) = log10(E/TeV) - bin centre\n")
232   - fprf.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n")
233   - fprf.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n")
234   - fprf.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n")
235   - fprf.write(" 5) Fractional Energy Resolution (rms)\n")
236   - fprf.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n")
237   - fprf.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n")
238   - fprf.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n")
239   - fprf.write(" systematics and statistics and at least 10 photons.\n")
  228 + # Write trailer
  229 + fprf.write("---------------------------------------------\n")
  230 + fprf.write("Notes\n")
  231 + fprf.write(" 1) log(E) = log10(E/TeV) - bin centre\n")
  232 + fprf.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n")
  233 + fprf.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n")
  234 + fprf.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n")
  235 + fprf.write(" 5) Fractional Energy Resolution (rms)\n")
  236 + fprf.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n")
  237 + fprf.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n")
  238 + fprf.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n")
  239 + fprf.write(" systematics and statistics and at least 10 photons.\n")
240 240  
241   - # Close files
242   - fprf.close()
243   - fbkg.close()
  241 + # Close files
  242 + fprf.close()
  243 + fbkg.close()
244 244  
245   - # Return
246   - return
  245 + # Return
  246 + return
247 247  
248 248  
249 249 #==========================#
250 250 # Main routine entry point #
251 251 #==========================#
252 252 if __name__ == '__main__':
253   - """
254   - """
255   - # Get ROOT filename from command line argument
256   - if (len(sys.argv) > 1):
257   - filename = sys.argv[1]
258   - else:
259   - #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEPerformanceBCE_Nov2010"
260   - #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_subarrays_root"
261   - caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_50deg_root"
  253 + """
  254 + """
  255 + # Get ROOT filename from command line argument
  256 + if (len(sys.argv) > 1):
  257 + filename = sys.argv[1]
  258 + else:
  259 + #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEPerformanceBCE_Nov2010"
  260 + #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_subarrays_root"
  261 + caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_50deg_root"
262 262  
263   - # Get list of all performance files
264   - files = glob.glob(caldb+"/*.root")
265   - files.sort()
266   -
267   - # Loop over all files
268   - for file in files:
269   -
270   - # Get filename
271   - head, tail = os.path.split(file)
272   -
273   - # Build name
274   - name = tail.strip(".root")
275   - if name.find("kb_") != -1:
276   - name = name.replace("_20deg", "")
277   - elif name.find("IFAE_") != -1:
278   - name = name.replace("Subarray", "")
279   - name = name.replace("IFAE_", "")
280   - name = "IFAE_"+name
281   - i = name.find("hours")
282   - if i != -1:
283   - name = name[0:i+1]
284   -
285   - # Build output filenames
286   - prfname = name + ".dat"
287   - bkgname = "bkg_"+name+".txt"
288   - #print prfname, bkgname
  263 + # Get list of all performance files
  264 + files = glob.glob(caldb+"/*.root")
  265 + files.sort()
  266 +
  267 + # Loop over all files
  268 + for file in files:
  269 +
  270 + # Get filename
  271 + head, tail = os.path.split(file)
  272 +
  273 + # Build name
  274 + name = tail.strip(".root")
  275 + if name.find("kb_") != -1:
  276 + name = name.replace("_20deg", "")
  277 + elif name.find("IFAE_") != -1:
  278 + name = name.replace("Subarray", "")
  279 + name = name.replace("IFAE_", "")
  280 + name = "IFAE_"+name
  281 + i = name.find("hours")
  282 + if i != -1:
  283 + name = name[0:i+1]
  284 +
  285 + # Build output filenames
  286 + prfname = name + ".dat"
  287 + bkgname = "bkg_"+name+".txt"
  288 + #print prfname, bkgname
289 289  
290   - # Create files
291   - make_files(file, prfname, bkgname)
  290 + # Create files
  291 + make_files(file, prfname, bkgname)
... ...
src/ctlike/ctlike.cpp
  View file @89ab90b
... ... @@ -356,15 +356,18 @@ void ctlike::get_parameters(void)
356 356 // Get event file name
357 357 std::string filename = (*this)["infile"].filename();
358 358  
359   - // Try first to open as unbinned FITS file
  359 + // Try first to open as FITS file
360 360 try {
361 361  
362 362 // Determine whether FITS file has EVENTS extension
  363 + /*
363 364 GFits fits(filename);
364 365 bool is_unbinned = fits.contains("EVENTS");
365 366 fits.close();
  367 + */
366 368  
367 369 // If FITS file has an events header then load as unbinned
  370 + /*
368 371 if (is_unbinned) {
369 372 obs.load_unbinned(filename);
370 373 }
... ... @@ -373,8 +376,12 @@ void ctlike::get_parameters(void)
373 376 else {
374 377 obs.load_binned(filename);
375 378 }
  379 + */
376 380  
377   - // Get task parameters for event file loading
  381 + // Load data
  382 + obs.load(filename);
  383 +
  384 + // Get other task parameters
378 385 m_stat = gammalib::toupper((*this)["stat"].string());
379 386 m_caldb = (*this)["caldb"].string();
380 387 m_irf = (*this)["irf"].string();
... ... @@ -408,7 +415,45 @@ void ctlike::get_parameters(void)
408 415 // Load observations from XML file
409 416 m_obs.load(filename);
410 417  
411   - }
  418 + // Check if all observations have response information. If
  419 + // not, get the calibration database parameters and set
  420 + // the response properly
  421 + bool asked_for_response = false;
  422 + for (int i = 0; i < m_obs.size(); ++i) {
  423 +
  424 + // Get CTA observation
  425 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  426 +
  427 + // Continue only if observation is a CTA observation
  428 + if (obs != NULL) {
  429 +
  430 + // If there is no effective area information then
  431 + // set response from task parameters
  432 + if (obs->response().aeff() == NULL) {
  433 +
  434 + // Get calibration parameters
  435 + if (!asked_for_response) {
  436 + m_caldb = (*this)["caldb"].string();
  437 + m_irf = (*this)["irf"].string();
  438 + asked_for_response = true;
  439 + }
  440 +
  441 + // Set response
  442 + GCaldb caldb;
  443 + if (gammalib::dir_exists(m_caldb)) {
  444 + caldb.rootdir(m_caldb);
  445 + }
  446 + else {
  447 + caldb.open("cta", m_caldb);
  448 + }
  449 + obs->response(m_irf, caldb);
  450 +
  451 + } // endif: there was no response
  452 + } // endif: observation was a CTA observation
  453 +
  454 + } // endfor: looped over observations
  455 +
  456 + } // endcatch: file was an XML file
412 457  
413 458 } // endif: there was no observation in the container
414 459  
... ...
src/ctlike/ctlike.par
  View file @89ab90b
... ... @@ -17,14 +17,13 @@
17 17 #
18 18 # General parameters
19 19 #===================
20   -infile, f, a, "events.fits",,,"Event list, counts map or observation definition file"
  20 +infile, f, a, "events.fits",,, "Event list, counts map or observation definition file"
21 21 stat, s, h, "POISSON",,, "Optimization statistics"
22 22 refit, b, h, no,,, "Do refitting?"
23   -edisp, b, h, no,,, "Apply energy dispersion?"
24   -caldb, s, a, "$CTOOLS/share/caldb/data/cta/dummy",,, "Calibration database"
  23 +caldb, s, a, "dummy",,, "Calibration database"
25 24 irf, s, a, "cta_dummy_irf",,, "Instrument response function"
26 25 srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
27   -outmdl, f, a, "crab_results.xml",,, "Source model output file"
  26 +outmdl, f, a, "crab_results.xml",,, "Source model output file"
28 27 #
29 28 # Standard parameters
30 29 #====================
... ...
src/ctobssim/ctobssim.par
  View file @89ab90b
... ... @@ -18,12 +18,12 @@
18 18 # File information
19 19 #=================
20 20 infile, f, a, "$CTOOLS/share/models/crab.xml",,, "Model"
21   -outfile,f, a, "events.fits",,, "Output event data file or observation definition file"
22   -prefix, s, h, "sim_events_",,, "Prefix for event lists in observation definition file"
  21 +outfile,f, a, "events.fits",,, "Output event data file or observation definition file"
  22 +prefix, s, h, "sim_events_",,, "Prefix for event lists in observation definition file"
23 23 #
24 24 # Observation parameters
25 25 #=======================
26   -caldb, s, a, "$CTOOLS/share/caldb/data/cta/dummy",,, "Calibration database"
  26 +caldb, s, a, "dummy",,, "Calibration database"
27 27 irf, s, a, "cta_dummy_irf",,, "Instrument response function"
28 28 #
29 29 # Simulation parameters
... ... @@ -31,7 +31,7 @@ irf, s, a, &quot;cta_dummy_irf&quot;,,, &quot;Instrument response function&quot;
31 31 seed, i, h, 1, , , "Random number generator seed"
32 32 ra, r, a, 83.63, 0,360, "RA of pointing (degrees)"
33 33 dec, r, a, 22.01, -90, 90, "Dec of pointing (degrees)"
34   -rad, r, a, 10.0, 0,180, "Radius of FOV (degrees)"
  34 +rad, r, a, 5.0, 0,180, "Radius of FOV (degrees)"
35 35 tmin, r, a, 0.0, 0, , "Start time (MET in s)"
36 36 tmax, r, a, 1800.0, 0, , "End time (MET in s)"
37 37 emin, r, a, 0.1, 0, , "Lower energy limit (TeV)"
... ...