diff --git a/ChangeLog b/ChangeLog index 67d9a53..c637531 100644 --- a/ChangeLog +++ b/ChangeLog @@ -6,12 +6,13 @@ Support energy dispersion in ctobssim, ctlike and ctmodel -2014-02-10 Juergen Knoedlseder +2014-02-20 Juergen Knoedlseder * Bug fix Version 0.7.2 released ================================ Add logL and event statistics to cspull output file + Properly close temporary file in ctselect (#1149) 2014-02-05 Juergen Knoedlseder diff --git a/NEWS b/NEWS index a300548..ce93dfa 100644 --- a/NEWS +++ b/NEWS @@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-08-00 Author(s) : Juergen Knoedlseder -10 February 2014 +20 February 2014 1. Introduction @@ -31,7 +31,8 @@ dispersion can be switched on or off using the hidden edisp parameter. 5. ctselect - CTA event selection --------------------------------- -None +The temporary file is now properly close to avoid file opening errors +when iterating over ctselect (#1149). 6. ctbin - CTA event binning diff --git a/caldb/Makefile.am b/caldb/Makefile.am index 0dd7d27..4d26926 100644 --- a/caldb/Makefile.am +++ b/caldb/Makefile.am @@ -1,10 +1,12 @@ ## Process this file with automake to produce Makefile.in # Calibration files to distribute -irfs = cta_dummy_irf.dat +caldb_indx = caldb.indx +caldb_irf = cta_dummy_irf.dat -caldbdir = $(datadir)/caldb/data/cta/dummy -caldb_DATA = $(irfs) - -# Extra files to be distributed in the tar file -EXTRA_DIST = $(irfs) +# Distribute calibration database files +EXTRA_DIST = $(caldb_indx) $(caldb_irf) +caldbindxdir = $(datadir)/caldb/data/cta/dummy +caldbindx_DATA = $(caldb_indx) +caldbirfdir = $(datadir)/caldb/data/cta/dummy/bcf +caldbirf_DATA = $(caldb_irf) diff --git a/caldb/caldb.indx b/caldb/caldb.indx new file mode 100644 index 0000000..5aa32fc --- /dev/null +++ b/caldb/caldb.indx @@ -0,0 +1 @@ +SIMPLE = T / file does conform to FITS standard BITPIX = 8 / number of bits per data pixel NAXIS = 0 / number of data axes EXTEND = T / FITS dataset may contain extensions COMMENT FITS (Flexible Image Transport System) format is defined in 'AstronomyCOMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H COMMENT FITS (Flexible Image Transport System) format is defined in 'AstronomyCOMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H END XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / 8-bit bytes NAXIS = 2 / 2-dimensional binary table NAXIS1 = 945 / width of table in bytes NAXIS2 = 4 / number of rows in table PCOUNT = 0 / size of special data area GCOUNT = 1 / one data group (required keyword) TFIELDS = 18 / number of fields in each row TTYPE1 = 'TELESCOP' / label for field 1 TFORM1 = '10A ' / data format of field: ASCII Character TTYPE2 = 'INSTRUME' / label for field 2 TFORM2 = '10A ' / data format of field: ASCII Character TTYPE3 = 'DETNAM ' / label for field 3 TFORM3 = '20A ' / data format of field: ASCII Character TTYPE4 = 'FILTER ' / label for field 4 TFORM4 = '10A ' / data format of field: ASCII Character TTYPE5 = 'CAL_DEV ' / label for field 5 TFORM5 = '20A ' / data format of field: ASCII Character TTYPE6 = 'CAL_DIR ' / label for field 6 TFORM6 = '70A ' / data format of field: ASCII Character TTYPE7 = 'CAL_FILE' / label for field 7 TFORM7 = '40A ' / data format of field: ASCII Character TTYPE8 = 'CAL_CLAS' / label for field 8 TFORM8 = '3A ' / data format of field: ASCII Character TTYPE9 = 'CAL_DTYP' / label for field 9 TFORM9 = '4A ' / data format of field: ASCII Character TTYPE10 = 'CAL_CNAM' / label for field 10 TFORM10 = '20A ' / data format of field: ASCII Character TTYPE11 = 'CAL_CBD ' / label for field 11 TFORM11 = '630A70 ' / data format of field: ASCII Character TTYPE12 = 'CAL_XNO ' / label for field 12 TFORM12 = '1I ' / data format of field: 2-byte INTEGER TTYPE13 = 'CAL_VSD ' / label for field 13 TFORM13 = '10A ' / data format of field: ASCII Character TTYPE14 = 'CAL_VST ' / label for field 14 TFORM14 = '8A ' / data format of field: ASCII Character TTYPE15 = 'REF_TIME' / label for field 15 TFORM15 = '1D ' / data format of field: 8-byte DOUBLE TTYPE16 = 'CAL_QUAL' / label for field 16 TFORM16 = '1I ' / data format of field: 2-byte INTEGER TTYPE17 = 'CAL_DATE' / label for field 17 TFORM17 = '8A ' / data format of field: ASCII Character TTYPE18 = 'CAL_DESC' / label for field 18 TFORM18 = '70A ' / data format of field: ASCII Character EXTNAME = 'CIF ' / name of this extension CIFVERSN= '1992a ' / Version of CIF format HISTORY File modified by user 'jurgen' with fv on 2014-02-20T15:10:40 END CTA DUMMY NONE NONE ONLINE data/cta/dummy/bcf cta_dummy_irf.dat BCFDATAEFF_AREA NAME(cta_dummy_irf) VERSION(Prod1) CLASS(BEST) ANALYSIS(MPIK) NONE NONE NONE NONE NONE 2014-01-0100:00:00@é+14/01/01CTA effective area CTA DUMMY NONE NONE ONLINE data/cta/dummy/bcf cta_dummy_irf.dat BCFDATARPSF NAME(cta_dummy_irf) VERSION(Prod1) CLASS(BEST) ANALYSIS(MPIK) NONE NONE NONE NONE NONE 2014-01-0100:00:00@é+14/01/01CTA point spread function CTA DUMMY NONE NONE ONLINE data/cta/dummy/bcf cta_dummy_irf.dat BCFDATAEDISP NAME(cta_dummy_irf) VERSION(Prod1) CLASS(BEST) ANALYSIS(MPIK) NONE NONE NONE NONE NONE 2014-01-0100:00:00@é+14/01/01CTA energy dispersion CTA DUMMY NONE NONE ONLINE data/cta/dummy/bcf cta_dummy_irf.dat BCFDATABGD NAME(cta_dummy_irf) VERSION(Prod1) CLASS(BEST) ANALYSIS(MPIK) NONE NONE NONE NONE NONE 2014-01-0100:00:00@é+14/01/01CTA background \ No newline at end of file diff --git a/doc/source/_templates/index.html b/doc/source/_templates/index.html index ff701d0..2597495 100644 --- a/doc/source/_templates/index.html +++ b/doc/source/_templates/index.html @@ -34,7 +34,8 @@ CTA format), Fermi/LAT and COMPTEL data. {%endtrans%}

  • {%trans%}Installation{%endtrans%}
    • -
  • {%trans%}Quickstart{%endtrans%}
    • +
  • {%trans%}Quickstart{%endtrans%}
    • +
  • {%trans%}CTA response functions{%endtrans%}
  • {%trans%}User Manual{%endtrans%}
  • {%trans%}ctools reference{%endtrans%}
  • {%trans%}Source code documentation (Doxygen){%endtrans%}
    • diff --git a/doc/source/conf.py b/doc/source/conf.py index c7c62bd..0aa0c60 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -51,7 +51,7 @@ copyright = u'2014, ctools developers' version = '0.7' # The full version, including alpha/beta/rc tags. -release = '0.7.1' +release = '0.7.2' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. @@ -165,7 +165,7 @@ html_show_copyright = False #html_file_suffix = None # Output file base name for HTML help builder. -htmlhelp_basename = 'GammaLibdoc' +htmlhelp_basename = 'ctoolsdoc' # -- Options for LaTeX output -------------------------------------------------- diff --git a/doc/source/user_manual/cntmap-fits.jpg b/doc/source/user_manual/getting_started/cntmap-fits.jpg index f9c73aa..f9c73aa 100644 --- a/doc/source/user_manual/cntmap-fits.jpg +++ b/doc/source/user_manual/getting_started/cntmap-fits.jpg diff --git a/doc/source/user_manual/cntmap-map.jpg b/doc/source/user_manual/getting_started/cntmap-map.jpg index b22de89..b22de89 100644 --- a/doc/source/user_manual/cntmap-map.jpg +++ b/doc/source/user_manual/getting_started/cntmap-map.jpg diff --git a/doc/source/user_manual/dskeys.jpg b/doc/source/user_manual/getting_started/dskeys.jpg index 0185466..0185466 100644 --- a/doc/source/user_manual/dskeys.jpg +++ b/doc/source/user_manual/getting_started/dskeys.jpg diff --git a/doc/source/user_manual/getting_started/faq.rst b/doc/source/user_manual/getting_started/faq.rst new file mode 100644 index 0000000..d6444b9 --- /dev/null +++ b/doc/source/user_manual/getting_started/faq.rst @@ -0,0 +1,4 @@ +Frequently asked questions +-------------------------- + +To be written ... diff --git a/doc/source/user_manual/getting_started/index.rst b/doc/source/user_manual/getting_started/index.rst new file mode 100644 index 0000000..d5ddba6 --- /dev/null +++ b/doc/source/user_manual/getting_started/index.rst @@ -0,0 +1,12 @@ +Getting started +=============== + +.. toctree:: + :maxdepth: 2 + + quickstart + response + models + python + faq + diff --git a/doc/source/user_manual/getting_started/irfs.jpg b/doc/source/user_manual/getting_started/irfs.jpg new file mode 100644 index 0000000..c60bb63 --- /dev/null +++ b/doc/source/user_manual/getting_started/irfs.jpg diff --git a/doc/source/user_manual/getting_started/models.rst b/doc/source/user_manual/getting_started/models.rst new file mode 100644 index 0000000..d9fa054 --- /dev/null +++ b/doc/source/user_manual/getting_started/models.rst @@ -0,0 +1,5 @@ +Models available in ctools +-------------------------- + +To be written ... (you may read +http://gammalib.sourceforge.net/user_manual/modules/model.html#spectral-models for the time being). diff --git a/doc/source/user_manual/getting_started/python.rst b/doc/source/user_manual/getting_started/python.rst new file mode 100644 index 0000000..2c75456 --- /dev/null +++ b/doc/source/user_manual/getting_started/python.rst @@ -0,0 +1,4 @@ +Using ctools from Python +------------------------ + +To be written ... diff --git a/doc/source/user_manual/quickstart.rst b/doc/source/user_manual/getting_started/quickstart.rst index 875416e..13b5a5f 100644 --- a/doc/source/user_manual/quickstart.rst +++ b/doc/source/user_manual/getting_started/quickstart.rst @@ -1,6 +1,3 @@ -Getting started -=============== - A quick ctools tutorial ----------------------- @@ -641,22 +638,3 @@ the same events for both analyses. 2014-01-09T21:47:54: +==============+ 2014-01-09T21:47:54: 2014-01-09T21:47:54: Application "ctlike" terminated after 6 wall clock seconds, consuming 0.163384 seconds of CPU time. - - -Models available in ctools --------------------------- - -To be written ... (you may read -http://gammalib.sourceforge.net/user_manual/modules/model.html#spectral-models for the time being). - - -Using ctools from Python ------------------------- - -To be written ... - - -Frequently asked questions --------------------------- - -To be written ... diff --git a/doc/source/user_manual/getting_started/response.rst b/doc/source/user_manual/getting_started/response.rst new file mode 100644 index 0000000..e65d5c6 --- /dev/null +++ b/doc/source/user_manual/getting_started/response.rst @@ -0,0 +1,395 @@ +CTA Instrument Response Functions +--------------------------------- + +.. note :: + + Instrument response functions are formally not part of ctools as they + should be provided by the instrument teams. The GammaLib framework on + which ctools are built comes with a set of basic response functions, but + for getting the latest instrument response functions you should contact + the relevant instrument teams. CTA consortium members should + `download + `_ + the latest calibration database containing Prod1 and Prod2 instrument + response functions from the consortium Sharepoint site (requires + password). + + +What are instrument response functions? +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The instrument response functions provide a mathematical description that +links the measured quantities :math:`\vec{d}` of an event to the physical +quantities :math:`\vec{p}` of the incident photon. The following figure +illustrates this relationship: + +.. figure:: irfs.jpg + :width: 70% + +:math:`I(\vec{p})` is the gamma-ray intensity arriving at Earth as a +function of photon properties :math:`\vec{p}` +(which usually are true photon energy, true photon incident direction, +and true photon arrival time), +while :math:`e(\vec{d})` is the expected event rate as function of event +properties :math:`\vec{d}` (which usually are the measured photon energy, +measured or reconstructued photon incident direction, and measured photon +arrival time). The expected event rate is obtained by integrating the +product of +the instrumental response function :math:`R(\vec{d}|\vec{p},\vec{a}`) +and the emitted intensity :math:`I(\vec{p})` over the photon properties +:math:`\vec{p}`. +The argument :math:`\vec{a}` in the response function comprises any +auxiliary parameter on which the response function may depend on (e.g. +pointing direction, triggered telescopes, optical efficiencies, +atmospheric conditions, etc.). All these quantities and hence the +instrument response function may depend on time. + + +Data formats +~~~~~~~~~~~~ + +So far the data format of the CTA instrument response functions has not +yet been settled, and information about instrument response functions is +provided in various ways. Three types of data formats are supported: + +- :ref:`sec_cta_rsp_perftable` + +- :ref:`sec_cta_rsp_xspec` + +- :ref:`sec_cta_rsp_rsptable` + +All formats provide the instrument response as function of true (and +sometimes also measured) photon energy, typically from about 20 GeV to +about 125 TeV. +Performance tables provide the instrument response for on-axis sources. +ARF, RMF and PSF files provide the instrument response for a specific +source position (and region) within the field of view. +Response tables provide the instrument response as function of position in +the field of view. Response tables are thus the most universal form of +instrument response functions available, and we recommend using this form +for any analysis. + +.. note :: + + Instrument response functions are so far only provided for a fixed + zenith angle of 20 deg. Therefore, no zenith or azimuth angle dependence has + been implemented so far. + + +Installing the CTA calibration database +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +After `downloading +`_ +the latest calibration database (only possible for CTA consortium members), +the database is installed using + +.. code-block:: bash + + $ [sudo] tar -C $CALDB -zxvf cta-caldb-20140220.tar.gz + +or + +.. code-block:: bash + + $ [sudo] tar -C $CALDB -xvf cta-caldb-20140220.tar + +depending on whether the database has been retrieved as a gzipped file or +not. Depending on your access rights to the ``$CALDB`` directory, +installation of the calibration database may require root privileges +(type ``sudo`` in this case, otherwise omit this part of the command). + +The calibration database must be installed into the directory to which you +``CALDB`` environment variable points. By default, ctools sets this +enviroment variable to + +.. code-block:: bash + + $CTOOLS/share/caldb + +but it may be that due to some other analysis software on your system the +``CALDB`` environment variable points to another directory. If in doubt, type + +.. code-block:: bash + + echo $CALDB + +to find out to which directory your ``CALDB`` environment variable points. +You may always overwrite this setting using + +.. code-block:: bash + + export CALDB=/my/preferred/caldb/directory + +which you can add to your ``.bashrc`` file to set the directory +permanently. +ctools will use the ``CALDB`` environment variable to find out where the +calibration database is located. If needed, this mechanism can be +circumvented (see :ref:`sec_cta_rsp_abspath`). + + +Specifying the CTA Instrument Response Functions +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The specification of the CTA Instrument Response Functions depends on the +way how ctools are used. Common to all methods is that the IRFs are +defined by a response name and a calibration database name. + +There are different means to specify the CTA Instrument Response Functions +when using ctools, and the following section describe the + + +Specifying the response function as input parameters +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +ctools that require instrument response functions have two parameters +to specify the calibration database name and the response function name. +The following example shows a ``ctobssim`` run using the ``dummy`` +calibration database and the ``cta_dummy_irf`` response function: + +.. code-block:: bash + + $ ctobssim + Model [$CTOOLS/share/models/crab.xml] + Calibration database [dummy] + Instrument response function [cta_dummy_irf] + RA of pointing (degrees) (0-360) [83.63] + Dec of pointing (degrees) (-90-90) [22.01] + Radius of FOV (degrees) (0-180) [5.0] + Start time (MET in s) (0) [0.0] + End time (MET in s) (0) [1800.0] + Lower energy limit (TeV) (0) [0.1] + Upper energy limit (TeV) (0) [100.0] + Output event data file or observation definition file [events.fits] + +Running the other tools is equivalent. + + +Specifying the response function in an observation definition file +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In the above example, only a single global response function can be +used for all CTA observations. If you need to specify response functions +per observation you can add the information directly in the XML observation +definition file: + +.. code-block:: xml + + + + + + + + +The ``Calibration`` parameter specifies the calibration database and +response name. You can then pass this file directly to, e.g., ``ctlike``: + +.. code-block:: bash + + $ ctlike + Event list, counts map or observation definition file [events.fits] obs_rsp.xml + Source model [$CTOOLS/share/models/crab.xml] + Source model output file [crab_results.xml] + +Note that ``ctlike`` does not ask for the calibration database and +response name as it found the relevant information in the XML file. + +.. _sec_cta_rsp_abspath: + + +Specifying individual instrument response files +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +If you need even more control over individual response files, you can +specify them individually in the XML observation file as follows: + +.. code-block:: xml + + + + + + + + + + + + +From within a Python script +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The following example illustrates how to set the calibration database +and response name from within Python: + +.. code-block:: python + + import gammalib + obs = gammalib.GCTAObservation() + caldb = gammalib.GCaldb("cta", "dummy") + irf = "cta_dummy_irf" + obs.response(irf, caldb) + +The calibration database is set by creating a ``GCaldb`` object. The +constructor takes as argument the mission (always ``cta``) and the +database name, in our case ``dummy``. The response function is then set +by passing the response name (here ``cta_dummy_irf``) and the calibration +database object to the ``response`` method. + +Alternatively, you can specify a specific response file by using + +.. code-block:: python + + import gammalib + obs = gammalib.GCTAObservation() + caldb = gammalib.GCaldb("$CALDB/data/cta/dummy/bcf") + irf = "cta_dummy_irf.dat" + obs.response(irf, caldb) + +Here, the directory in which the calibration file resides is specified +as argument to the ``GCaldb`` constructor, and the filename is passed +as response function name to the ``response`` method. + + +What response functions are available? +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +To see what response functions are available on your system, type + +.. code-block:: bash + + $ cscaldb + +This creates the file ``cscaldb.log`` in your current working directory +that contains information about the available response functions. The +format of the file is: + +.. code-block:: bash + + 2014-02-20T17:15:21: +==============+ + 2014-02-20T17:15:21: | Mission: cta | + 2014-02-20T17:15:21: +==============+ + 2014-02-20T17:15:21: === Calibration: aar === + 2014-02-20T17:15:21: DESY20140105_50h + 2014-02-20T17:15:21: DESY20140105_50h_0deg + 2014-02-20T17:15:21: DESY20140105_50h_180deg + 2014-02-20T17:15:21: + 2014-02-20T17:15:21: === Calibration: aar500 === + 2014-02-20T17:15:21: DESY20140105_50h + 2014-02-20T17:15:21: DESY20140105_50h_0deg + 2014-02-20T17:15:21: DESY20140105_50h_180deg + +The name after ``Calibration`` is what we called earlier the calibration +database, and the names that are following each +``=== Calibration: xxx ===`` header are the response file names. + +To get the same information also onto the console, type + +.. code-block:: bash + + $ cscaldb debug=yes + + +Data format details +~~~~~~~~~~~~~~~~~~~ + +.. _sec_cta_rsp_perftable: + +Performance tables +^^^^^^^^^^^^^^^^^^ + +In the early days, the instrument performances derived from Monte-Carlo +simulations have been summarised in what we call here Performance +Tables, which are ASCII files that contain as function of energy the +on-axis performance parameters of CTA, such as effective area, point spread +function containment radius, energy resolution, background count rate and +differential sensitivity. + +Below an example of a CTA performance table:: + + log(E) Area r68 r80 ERes. BG Rate Diff Sens + -1.7 261.6 0.3621 0.4908 0.5134 0.0189924 6.88237e-11 + -1.5 5458.2 0.2712 0.3685 0.4129 0.1009715 1.72717e-11 + -1.3 15590.0 0.1662 0.2103 0.2721 0.0575623 6.16963e-12 + -1.1 26554.1 0.1253 0.1567 0.2611 0.0213008 2.89932e-12 + -0.9 52100.5 0.1048 0.1305 0.1987 0.0088729 1.39764e-12 + -0.7 66132.1 0.0827 0.1024 0.1698 0.0010976 6.03531e-13 + -0.5 108656.8 0.0703 0.0867 0.1506 0.0004843 3.98147e-13 + -0.3 129833.0 0.0585 0.0722 0.1338 0.0001575 3.23090e-13 + -0.1 284604.3 0.0531 0.0656 0.1008 0.0001367 2.20178e-13 + 0.1 263175.3 0.0410 0.0506 0.0831 0.0000210 1.87452e-13 + 0.3 778048.6 0.0470 0.0591 0.0842 0.0000692 1.53976e-13 + 0.5 929818.8 0.0391 0.0492 0.0650 0.0000146 1.18947e-13 + 0.7 1078450.0 0.0335 0.0415 0.0541 0.0000116 1.51927e-13 + 0.9 1448579.1 0.0317 0.0397 0.0516 0.0000047 1.42439e-13 + 1.1 1899905.0 0.0290 0.0372 0.0501 0.0000081 1.96670e-13 + 1.3 2476403.8 0.0285 0.0367 0.0538 0.0000059 2.20695e-13 + 1.5 2832570.6 0.0284 0.0372 0.0636 0.0000073 3.22523e-13 + 1.7 3534065.3 0.0290 0.0386 0.0731 0.0000135 4.84153e-13 + 1.9 3250103.4 0.0238 0.0308 0.0729 0.0000044 6.26265e-13 + 2.1 3916071.6 0.0260 0.0354 0.0908 0.0000023 7.69921e-13 + --------------------------------------------- + 1) log(E) = log10(E/TeV) - bin centre + 2) Eff Area - in square metres after background cut (no theta cut) + 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees + 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees + 5) Fractional Energy Resolution (rms) + 6) BG Rate - inside point-source selection region - post call cuts - in Hz + 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE + - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including + systematics and statistics and at least 10 photons. + + +.. _sec_cta_rsp_xspec: + +ARF, RMF and PSF files +^^^^^^^^^^^^^^^^^^^^^^ + +Instrument response information for the first CTA Data Challenge (1DC) has been +provided in a format that was heavily inspired by the +`RMF and ARF file formats +`_ +that have been introduced by +`HEASARC `_ +for spectral analysis of X-ray data. + +The Redistribution Matrix File (RMF) describes how an incoming photon with +a given true energy is redistributed in measured energy. In other words, it +describes the energy dispersion of the instrument. The RMF is organised as +a two-dimensional matrix, with the first axis being in true energies while +the second axis represents the measured energies. + +The Ancillary Response File (ARF) describes the sensitivity of the +instrument to photons of a given true energy. The ARF gives the effective +area of the detector system after applying an event selection cut +(theta cut). For computing the ARF, some knowledge of the Point Spread +Function (PSF) is needed, so that the fraction of photons that falls into +the event selection region for a given source can be properly estimated. +Once this is done, no further PSF information is needed for the analysis. + +ctools however relies on the full modelling of the instrument, including +the Point Spread Function, and hence, PSF information has also been +provided for 1DC. Two data formats have been used for this: a first that +is based on a simple one-dimensional vector, providing the width of a +2-dimensional Gaussian as function of energy; and a second that is based +on a three-component 2D Gaussian as function of energy and offset angle +in the camera system. While the former was implement by a simply FITS +table with two columns, the latter was implemented by +:ref:`sec_cta_rsp_rsptable`. + + +.. _sec_cta_rsp_rsptable: + +Response tables +^^^^^^^^^^^^^^^ + +The CTA response table class ``GCTAResponseTable`` provides a generic +handle for multi-dimensional response information. It is based on the +response format used for storing response information for the +*Fermi*/LAT telescope. In this format, all information is stored in +a single row of a FITS binary table. Each element of the row contains +a vector column, that describes the axes of the multi-dimensional response +cube and the response information. Note that this class may in the future +be promoted to the GammaLib core, as a similar class has been implemented +in the *Fermi*/LAT interface. diff --git a/doc/source/user_manual/index.rst b/doc/source/user_manual/index.rst index 079b194..0c5a40d 100644 --- a/doc/source/user_manual/index.rst +++ b/doc/source/user_manual/index.rst @@ -7,7 +7,7 @@ User Manual introduction getting - quickstart + getting_started/index enhancing reference expert diff --git a/scripts/Makefile.am b/scripts/Makefile.am index 810837b..c23790b 100644 --- a/scripts/Makefile.am +++ b/scripts/Makefile.am @@ -4,12 +4,12 @@ binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py # Scripts to be distributed in pythondir -pythonscripts = cspull.py cssens.py cstsdist.py \ +pythonscripts = cspull.py cssens.py cstsdist.py cscaldb.py \ obsutils.py \ lsmodel.py # Scripts for which symbolic links will be created in bindir -linkscripts = cspull cssens cstsdist \ +linkscripts = cspull cssens cstsdist cscaldb \ lsmodel # Distribute scripts diff --git a/scripts/cscaldb.py b/scripts/cscaldb.py new file mode 100755 index 0000000..dadc680 --- /dev/null +++ b/scripts/cscaldb.py @@ -0,0 +1,209 @@ +#! /usr/bin/env python +# ========================================================================== +# This script dumps all available calibrations into the console. +# +# Copyright (C) 2014 Juergen Knoedlseder +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import gammalib +import sys +import glob +import os + + +# ============ # +# cspull class # +# ============ # +class cscaldb(gammalib.GApplication): + """ + This class implements the pull distribution generation script. It derives + from the GammaLib::GApplication class which provides support for parameter + files, command line arguments, and logging. In that way the Python + script behaves just as a regular ctool. + """ + def __init__(self, *argv): + """ + Constructor. + """ + # Set name + self.name = "cscaldb" + self.version = "0.1.0" + + # Make sure that parfile exists + file = self.parfile() + + # Initialise application + if len(argv) == 0: + gammalib.GApplication.__init__(self, self.name, self.version) + elif len(argv) ==1: + gammalib.GApplication.__init__(self, self.name, self.version, *argv) + else: + raise TypeError("Invalid number of arguments given.") + + # Set logger properties + self.log_header() + self.log.date(True) + + # Return + return + + def __del__(self): + """ + Destructor. + """ + # Return + return + + def parfile(self): + """ + Check if parfile exists. If parfile does not exist then create a + default parfile. This kluge avoids shipping the cscript with a parfile. + """ + # Set parfile name + parfile = self.name+".par" + + try: + pars = gammalib.GApplicationPars(parfile) + except: + # Signal if parfile was not found + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") + + # Create default parfile + pars = gammalib.GApplicationPars() + pars.append_standard() + pars.save(parfile) + + # Return + return + + def get_parameters(self): + """ + Get parameters from parfile. + """ + # Return + return + + def execute(self): + """ + Execute the script. + """ + # Run the script + self.run() + + # Return + return + + def run(self): + """ + Run the script. + """ + # Switch screen logging on in debug mode + if self.logDebug(): + self.log.cout(True) + + # Get parameters + self.get_parameters() + + # Write input parameters into logger + if self.logTerse(): + self.log_parameters() + self.log("\n") + + # Get the calibration database path + caldb = gammalib.GCaldb() + + # Gather missions + missions = [] + paths = glob.glob(caldb.rootdir()+"/data/*") + paths.sort() + for path in paths: + head, tail = os.path.split(path) + missions.append(tail) + + # Loop over missions + for mission in missions: + + # Skip all non-CTA instruments + if mission != "cta": + continue + + # Write mission into logger + if self.logTerse(): + self.log("\n") + self.log.header1("Mission: "+mission) + + # Gather instruments + instruments = [] + paths = glob.glob(caldb.rootdir()+"/data/"+mission+"/*") + paths.sort() + for path in paths: + head, tail = os.path.split(path) + instruments.append(tail) + + # Loop over instruments + for instrument in instruments: + + # Write mission into logger + if self.logTerse(): + self.log.header3("Calibration: "+instrument) + + # Open calibration index file + cifname = caldb.rootdir()+"/data/"+mission+"/"+instrument+"/caldb.indx" + fits = gammalib.GFits(cifname) + cif = fits["CIF"] + cals = cif["CAL_CBD"] + + # Extract response names + names = [] + nrows = cals.length() + ncols = cals.number() + for row in range(nrows): + for col in range(ncols): + cal = cals.string(row, col) + istart = cal.find("NAME(") + if istart != -1: + istop = cal.find(")") + name = cal[5:istop] + if names.count(name) == 0: + names.append(name) + names.sort() + + # Print response name + if self.logTerse(): + for name in names: + self.log(name+"\n") + self.log("\n") + + # Return + return + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + """ + Dumps calibration database into log file. + """ + # Create instance of application + app = cscaldb(sys.argv) + + # Open logfile + app.logFileOpen() + + # Execute application + app.execute() + \ No newline at end of file diff --git a/scripts/cssens.py b/scripts/cssens.py index 242bee8..0bb8858 100755 --- a/scripts/cssens.py +++ b/scripts/cssens.py @@ -57,7 +57,7 @@ class cssens(GApplication): """ # Set name self.name = "cssens" - self.version = "0.3.0" + self.version = "0.4.0" # Make sure that parfile exists file = self.parfile() @@ -106,11 +106,11 @@ class cssens(GApplication): # Create default parfile pars = GApplicationPars() pars.append(GApplicationPar("outfile","f","h","sensitivity.dat","","","Output file name")) - pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database")) + pars.append(GApplicationPar("caldb","s","a","$CTOOLS/share/caldb/data/cta/dummy","","","Calibration database")) pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function")) pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)")) pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)")) - pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)")) + pars.append(GApplicationPar("bkg","s","a","$CTOOLS/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)")) pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)")) pars.append(GApplicationPar("rad","r","a","5.0","","","Radius of ROI (deg)")) pars.append(GApplicationPar("enumbins","i","h","0","","","Number of energy bins (0=unbinned)")) @@ -333,7 +333,7 @@ class cssens(GApplication): location.lb_deg(l, b) # Set source spectrum - spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV") + spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV")) if fitidx: spectrum["Index"].free() else: diff --git a/scripts/cta_kb_root2caldb.py b/scripts/cta_kb_root2caldb.py new file mode 100755 index 0000000..37c8244 --- /dev/null +++ b/scripts/cta_kb_root2caldb.py @@ -0,0 +1,545 @@ +#! /usr/bin/env python +# ========================================================================== +# This script generates CALDB compliant performance tables from MC ROOT +# files. +# +# Note that the effective area given in the ASCII file is the full effective +# area (i.e. the area before applying any theta cut). This is required by +# GammaLib as the maximum likelihood analysis does not imply any theta cut. +# +# Copyright (C) 2014 Juergen Knoedlseder +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +from ROOT import TFile, TH1F, TH2F +from datetime import datetime +import gammalib +import math +import sys +import glob +import os + + +# =========== # +# CalDB class # +# =========== # +class caldb(): + """ + """ + def __init__(self, inst, path=""): + """ + CALDB constructor. The constructor will create one CALDB entry for + an observation. + + Parameters: + inst - Instrument string. + Keywords: + path - CALDB Root path. + """ + # Store root path + self.path = path + + # Initialise some members + self.cif = None + self.irf = None + self.ea = None + self.psf = None + self.edisp = None + self.bgd = None + self.hdu_cif = None + + # Store UTC date of job execution + self.utc = datetime.utcnow().isoformat()[:19] + + # Set CALDB information + self.cal_tel = "CTA" + self.cal_inst = inst.upper() + self.cal_det = "NONE" + self.cal_flt = "NONE" + self.cal_class = "BCF" + self.cal_type = "DATA" + self.cal_qual = 0 # 0=good, 1=bad, 2=dubious, ... + self.cal_date = "14/01/01" + self.val_date = "2014-01-01" + self.val_time = "00:00:00" + self.ref_time = 51544.0 + self.cal_name = "NAME()" + self.cal_version = "VERSION(Prod1)" + self.cal_cut = "CLASS(BEST)" + self.cal_analysis = "ANALYSIS(MPIK)" + # + self.ea_name = "EFF_AREA" + self.ea_doc = "???" + self.ea_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] + self.ea_desc = "CTA effective area" + # + self.psf_name = "RPSF" + self.psf_doc = "???" + self.psf_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] + self.psf_desc = "CTA point spread function" + # + self.edisp_name = "EDISP" + self.edisp_doc = "???" + self.edisp_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] + self.edisp_desc = "CTA energy dispersion" + # + self.bgd_name = "BGD" + self.bgd_doc = "???" + self.bgd_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] + self.bgd_desc = "CTA background" + + # Create directory structure + self.make_dirs() + + # Return + return + + def __del__(self): + """ + Destructor. Close all FITS files. + """ + # Close calibration + self.close() + + # Return + return + + def make_dirs(self): + """ + Generate CALDB directory structure for one observation identifier. + The structure is given by + + data///bcf/ + + where is "cta" and is the instrument specified in the + CALDB constructor (the instrument may be used for different array + configurations). + + Parameters: + None + Keywords: + None + """ + # Set base path + self.base_dir = "data" + self.base_dir += "/"+self.cal_tel.lower() + self.base_dir += "/"+self.cal_inst.lower() + + # Set directory names + self.bcf_dir = self.base_dir+"/bcf" + + # Optionally prefix path + if len(self.path) > 0: + self.base_path = self.path+"/"+self.base_dir + self.bcf_path = self.path+"/"+self.bcf_dir + else: + self.base_path = self.base_dir + self.bcf_path = self.bcf_dir + + # Create directory structure + if not os.path.isdir(self.bcf_path): + os.makedirs(self.bcf_path) + + # Return + return + + def add(self, rspname, split=False, clobber=True): + """ + Add new calibration. The actual version will put + all calibrations in the same file, although each part of the response + function will have its own logical name. We can thus easily modify + the script to put each calibration information in a separate file. + + Parameters: + rspname - Response name + Keywords: + split - Split IRF over several files? + clobber - Overwrite existing files? + """ + # Set calibrate file names + if split: + self.ea_file = "ea_"+rspname+".dat" + self.psf_file = "psf_"+rspname+".dat" + self.edisp_file = "edisp_"+rspname+".dat" + self.bgd_file = "bgd_"+rspname+".dat" + else: + self.ea_file = "irf_"+rspname+".dat" + self.psf_file = "irf_"+rspname+".dat" + self.edisp_file = "irf_"+rspname+".dat" + self.bgd_file = "irf_"+rspname+".dat" + + # Open calibration database index + if self.cif == None: + self.cif = gammalib.GFits(self.base_path+"/caldb.indx", True) + + # If file has no CIF extension than create it now + try: + self.hdu_cif = self.cif.table("CIF") + except: + self.create_cif() + self.hdu_cif = self.cif.table("CIF") + + # Set response name + self.cal_name = "NAME("+rspname+")" + + # Return + return + + def close(self): + """ + Close calibration FITS files. + + Parameters: + None + Keywords: + None + """ + # Add information to CIF. We do this now as before we did not + # necessarily have all the information at hand (in particular + # about the boundaries) + #self.add_cif_info() + + # Close CIF + if self.cif != None: + self.cif.save(True) + self.cif.close() + self.cif = None + + # Free + self.cif = None + + # Return + return + + def create_cif(self): + """ + Create Calibration Index File (CIF) extension in CIF FITS file. + + Parameters: + None + Keywords: + None + """ + # Create binary table + table = gammalib.GFitsBinTable() + + # Set boundary + + # Attach columns. Reference: CAL/GEN/92-008 + table.append(gammalib.GFitsTableStringCol("TELESCOP", 0, 10)) + table.append(gammalib.GFitsTableStringCol("INSTRUME", 0, 10)) + table.append(gammalib.GFitsTableStringCol("DETNAM", 0, 20)) + table.append(gammalib.GFitsTableStringCol("FILTER", 0, 10)) + table.append(gammalib.GFitsTableStringCol("CAL_DEV", 0, 20)) + table.append(gammalib.GFitsTableStringCol("CAL_DIR", 0, 70)) + table.append(gammalib.GFitsTableStringCol("CAL_FILE", 0, 40)) + table.append(gammalib.GFitsTableStringCol("CAL_CLAS", 0, 3)) + table.append(gammalib.GFitsTableStringCol("CAL_DTYP", 0, 4)) + table.append(gammalib.GFitsTableStringCol("CAL_CNAM", 0, 20)) + table.append(gammalib.GFitsTableStringCol("CAL_CBD", 0, 70, 9)) + table.append(gammalib.GFitsTableShortCol("CAL_XNO", 0)) + table.append(gammalib.GFitsTableStringCol("CAL_VSD", 0, 10)) + table.append(gammalib.GFitsTableStringCol("CAL_VST", 0, 8)) + table.append(gammalib.GFitsTableDoubleCol("REF_TIME", 0)) + table.append(gammalib.GFitsTableShortCol("CAL_QUAL", 0)) + table.append(gammalib.GFitsTableStringCol("CAL_DATE", 0, 8)) + table.append(gammalib.GFitsTableStringCol("CAL_DESC", 0, 70)) + + # Set keywords. Reference: CAL/GEN/92-008 + table.extname("CIF") + table.card("CIFVERSN", "1992a", "Version of CIF format") + + # Attach table to FITS file. Note that at this moment the + # FITS table gets connected to the FITS file. Yet since + # nothing was yet written to the FITS file, we cannot read + # anything from it. + self.cif.append(table) + + # Return + return + + def add_cif_info(self): + """ + Add information to CIF extension. + + Parameters: + None + Keywords: + None + """ + # Set boundaries + self.ea_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] + self.psf_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] + self.edisp_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] + self.bgd_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis] + + # Append 4 rows to CIF extension + self.hdu_cif.append_rows(4) + + # Add generic information for these 4 rows + for i in range(4): + + # Set row number + row = i+self.hdu_cif.nrows()-4 + + # Set element + self.hdu_cif["TELESCOP"][row] = self.cal_tel + self.hdu_cif["INSTRUME"][row] = self.cal_inst + self.hdu_cif["DETNAM"][row] = self.cal_det + self.hdu_cif["FILTER"][row] = self.cal_flt + self.hdu_cif["CAL_DEV"][row] = "ONLINE" + self.hdu_cif["CAL_CLAS"][row] = self.cal_class + self.hdu_cif["CAL_DTYP"][row] = self.cal_type + self.hdu_cif["CAL_VSD"][row] = self.val_date + self.hdu_cif["CAL_VST"][row] = self.val_time + self.hdu_cif["REF_TIME"][row] = self.ref_time + self.hdu_cif["CAL_QUAL"][row] = self.cal_qual + self.hdu_cif["CAL_DATE"][row] = self.cal_date + + # Add effective area information + row = self.hdu_cif.nrows()-4 + self.hdu_cif["CAL_DIR"][row] = self.bcf_dir + self.hdu_cif["CAL_FILE"][row] = self.ea_file + self.hdu_cif["CAL_CNAM"][row] = self.ea_name + self.hdu_cif["CAL_DESC"][row] = self.ea_desc + self.hdu_cif["CAL_XNO"][row] = 1 + for n in range(9): + if n >= len(self.ea_bounds): + self.hdu_cif["CAL_CBD"][row, n] = "NONE" + else: + self.hdu_cif["CAL_CBD"][row, n] = self.ea_bounds[n] + + # Add point spread function information + row = self.hdu_cif.nrows()-3 + self.hdu_cif["CAL_DIR"][row] = self.bcf_dir + self.hdu_cif["CAL_FILE"][row] = self.psf_file + self.hdu_cif["CAL_CNAM"][row] = self.psf_name + self.hdu_cif["CAL_DESC"][row] = self.psf_desc + self.hdu_cif["CAL_XNO"][row] = 1 + for n in range(9): + if n >= len(self.psf_bounds): + self.hdu_cif["CAL_CBD"][row, n] = "NONE" + else: + self.hdu_cif["CAL_CBD"][row, n] = self.psf_bounds[n] + + # Add energy dispersion information + row = self.hdu_cif.nrows()-2 + self.hdu_cif["CAL_DIR"][row] = self.bcf_dir + self.hdu_cif["CAL_FILE"][row] = self.edisp_file + self.hdu_cif["CAL_CNAM"][row] = self.edisp_name + self.hdu_cif["CAL_DESC"][row] = self.edisp_desc + self.hdu_cif["CAL_XNO"][row] = 1 + for n in range(9): + if n >= len(self.edisp_bounds): + self.hdu_cif["CAL_CBD"][row, n] = "NONE" + else: + self.hdu_cif["CAL_CBD"][row, n] = self.edisp_bounds[n] + + # Add background information + row = self.hdu_cif.nrows()-1 + self.hdu_cif["CAL_DIR"][row] = self.bcf_dir + self.hdu_cif["CAL_FILE"][row] = self.bgd_file + self.hdu_cif["CAL_CNAM"][row] = self.bgd_name + self.hdu_cif["CAL_DESC"][row] = self.bgd_desc + self.hdu_cif["CAL_XNO"][row] = 1 + for n in range(9): + if n >= len(self.bgd_bounds): + self.hdu_cif["CAL_CBD"][row, n] = "NONE" + else: + self.hdu_cif["CAL_CBD"][row, n] = self.bgd_bounds[n] + + # Return + return + + def root2caldb(self, filename): + """ + Generate files. + + Parameters: + filename - Performance file name + """ + # Open ROOT performance file + rootfile = TFile(filename) + + # Set filenames + ea_filename = self.bcf_path+"/"+self.ea_file + psf_filename = self.bcf_path+"/"+self.psf_file + edisp_filename = self.bcf_path+"/"+self.edisp_file + bgd_filename = self.bcf_path+"/"+self.bgd_file + + # Open ASCII performance tables + if self.ea_file == self.psf_file and \ + self.ea_file == self.edisp_file and \ + self.ea_file == self.bgd_file : + split = False + fprf = open(ea_filename, "w") + files = [fprf] + else: + split = True + fea = open(ea_filename, "w") + fpsf = open(psf_filename, "w") + fedisp = open(edisp_filename, "w") + fbgd = open(bgd_filename, "w") + files = [fea, fpsf, fedisp, fbgd] + + # Write header + for file in files: + file.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n") + + # Get histograms. Konrad's root files do not have an EffectiveArea80 + # histogram, but his effective area is for 80% containment radius + # anyways, so we simply read this histogram into aeff80. + sens = TH1F() + bgrate = TH1F() + bgratesqdeg = TH1F() + aeff = TH1F() + aeff80 = TH1F() + angres = TH1F() + angres80 = TH1F() + eres = TH1F() + rootfile.GetObject("DiffSens", sens) + rootfile.GetObject("BGRate", bgrate) + rootfile.GetObject("BGRatePerSqDeg", bgratesqdeg) + rootfile.GetObject("EffectiveArea", aeff) + try: + rootfile.GetObject("EffectiveArea80", aeff80) + except: + rootfile.GetObject("EffectiveArea", aeff80) + rootfile.GetObject("AngRes", angres) + rootfile.GetObject("AngRes80", angres80) + rootfile.GetObject("ERes", eres) + + # Setup vectors + energies = bgrate.GetXaxis() + nbins_eng = energies.GetNbins() + for i in range(nbins_eng): + + # Get logE and background rate per squared degree + logE = energies.GetBinCenter(i+1) + bgd = bgrate.GetBinContent(i+1,1) + bgdsq = bgratesqdeg.GetBinContent(i+1,1) + area = aeff.GetBinContent(i+1,1) + area80 = aeff80.GetBinContent(i+1,1) + r68 = angres.GetBinContent(i+1,1) + r80 = angres80.GetBinContent(i+1,1) + dE = eres.GetBinContent(i+1,1) + diffs = sens.GetBinContent(i+1,1) + + # Skip any NaNs + if math.isnan(area80): + continue + + # Compute energy (in MeV) + energy = math.pow(10.0, logE)*1.0e6 + emin = math.pow(10.0, logE-0.1)*1.0e6 + emax = math.pow(10.0, logE+0.1)*1.0e6 + ewidth = emax-emin + + # Convert into background rate per steradian and MeV + bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth + + # Compute control background rate (this only works for KB files + # as they are for 80% containment radius). But we don't use this + # anyways, it's just for control on display ... + omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80))) + if omega > 0.0: + bkg_rate2 = bgd / omega / ewidth + else: + bkg_rate2 = 0.0 + + # Compute full effective area by dividing area80/0.80 + area = area80/0.80 + + # Write results in file + line = "%.1f %.1f %.4f %.4f %.4f %.5e %.5e" % \ + (logE, area, r68, r80, dE, bgd, diffs) + for file in files: + file.write(line+"\n") + #fbkg.write(str(energy)+" "+str(bkg_rate)+"\n") + + # Write trailer + for file in files: + file.write("---------------------------------------------\n") + file.write("Notes\n") + file.write(" 1) log(E) = log10(E/TeV) - bin centre\n") + file.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n") + file.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n") + file.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n") + file.write(" 5) Fractional Energy Resolution (rms)\n") + file.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n") + file.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n") + file.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n") + file.write(" systematics and statistics and at least 10 photons.\n") + + # Close files + for file in files: + file.close() + + # Close ROOT performance file + rootfile.Close() + + # Add information to CIF. We do this now as before we did not + # necessarily have all the information at hand (in particular + # about the boundaries) + self.add_cif_info() + + # Return + return + + +#==========================# +# Main routine entry point # +#==========================# +if __name__ == '__main__': + """ + """ + # Get ROOT filename from command line argument + if (len(sys.argv) > 1): + filename = sys.argv[1] + else: + db = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_20deg_root" + #db = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_50deg_root" + + # Get list of all performance files + files = glob.glob(db+"/*.root") + files.sort() + + # Allocate caldb + irf = caldb("kb") + + # Loop over all files + for file in files: + + # Get filename + head, tail = os.path.split(file) + + # Build name + name = tail.strip(".root") + if name.find("kb_") != -1: + name = name.replace("_20deg", "") + elif name.find("IFAE_") != -1: + name = name.replace("Subarray", "") + name = name.replace("IFAE_", "") + name = "IFAE_"+name + i = name.find("hours") + if i != -1: + name = name[0:i+1] + print(tail) + + # Open calibration file + irf.add(name) + + # Translate ROOT to CALDB information + irf.root2caldb(file) diff --git a/scripts/cta_root2caldb.py b/scripts/cta_root2caldb.py index 8618101..e406369 100755 --- a/scripts/cta_root2caldb.py +++ b/scripts/cta_root2caldb.py @@ -1319,9 +1319,9 @@ def set_prod1_ifae(): #=========================# # Set Prod2 DESY database # #=========================# -def set_prod2_desy(): +def set_prod2_desy_ignore(): """ - Set Prod2 DESY database. + Set Prod2 DESY database (ignore these files). """ # Set database attributes path_aar = "/project-data/cta/performance/prod2/Performance_DESY_20140131_Aar" @@ -1366,6 +1366,56 @@ def set_prod2_desy(): return db +#=========================# +# Set Prod2 DESY database # +#=========================# +def set_prod2_desy(): + """ + Set Prod2 DESY database. + """ + # Set database attributes + path_aar = "/project-data/cta/performance/prod2/Performance_DESY_20140131_Aar" + path_ten = "/project-data/cta/performance/prod2/Performance_DESY_20140131_Tenerife_50h" + rebin = False + psftype = "Gauss" + scale = 1.0 + + # Set database content + db = [{'inst': "aar", 'id': "DESY20140105_50h", + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale, + 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-NS.S.2a.180000s.root"}, + {'inst': "aar", 'id': "DESY20140105_50h_0deg", + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale, + 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-NS.S.2a.180000s.root"}, + {'inst': "aar", 'id': "DESY20140105_50h_180deg", + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale, + 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-NS.S.2a.180000s.root"}, + + {'inst': "aar500", 'id': "DESY20140105_50h", + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale, + 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-500m-NS.S.2a.180000s.root"}, + {'inst': "aar500", 'id': "DESY20140105_50h_0deg", + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale, + 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-NS.S.2a.180000s.root"}, + {'inst': "aar500", 'id': "DESY20140105_50h_180deg", + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale, + 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-NS.S.2a.180000s.root"}, + + {'inst': "tenerife", 'id': "DESY20140105_50h", + 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'scale': scale, + 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Tenerife-NS.N.2NN.180000s.root"}, + {'inst': "tenerife", 'id': "DESY20140105_50h_0deg", + 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'scale': scale, + 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"}, + {'inst': "tenerife", 'id': "DESY20140105_50h_180deg", + 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'scale': scale, + 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"} + ] + + # Return database + return db + + #==========================# # Main routine entry point # #==========================# @@ -1381,7 +1431,7 @@ if __name__ == '__main__': - EffectiveArea80_offaxis (2D) - Effective area for 80% source containment - AngRes_offaxis (2D) - Angular resolution 68% containment - AngRes80_offaxis (2D) - Angular resolution 80% containment - - BGRatePerSqDeg_offaxis (2D) - Background rate per square degree + - BGRatePerSqDeg_offaxis (2D) - Background rate per square degree [counts/s/deg2] - BGRate_offaxis (2D) - Background rate - EffectiveAreaEtrue_offaxis (2D) - Effective area in true energy (finer bins) - MigMatrix_offaxis (3D) diff --git a/scripts/cta_root2irf.py b/scripts/cta_root2irf.py index c12f7b4..466750e 100755 --- a/scripts/cta_root2irf.py +++ b/scripts/cta_root2irf.py @@ -38,58 +38,58 @@ import os # Compute Eq. (17) of Li & Ma (1983) # # ================================== # def sigma_lima(Non, Noff, alpha=0.2): - """ - Compute Eq. (17) of Li & Ma (1983). - - Parameters: - Non - Number of on counts - Noff - Number of off counts - Keywords: - alpha - Ratio of on-to-off exposure - """ - # Compute sensitivity - alpha1 = alpha + 1.0 - sum = Non + Noff - arg1 = Non / sum - arg2 = Noff / sum - term1 = Non * math.log((alpha1/alpha)*arg1) - term2 = Noff * math.log(alpha1*arg2) - sigma = math.sqrt(2.0 * (term1 + term2)) - - # Return sensitivity - return sigma + """ + Compute Eq. (17) of Li & Ma (1983). + + Parameters: + Non - Number of on counts + Noff - Number of off counts + Keywords: + alpha - Ratio of on-to-off exposure + """ + # Compute sensitivity + alpha1 = alpha + 1.0 + sum = Non + Noff + arg1 = Non / sum + arg2 = Noff / sum + term1 = Non * math.log((alpha1/alpha)*arg1) + term2 = Noff * math.log(alpha1*arg2) + sigma = math.sqrt(2.0 * (term1 + term2)) + + # Return sensitivity + return sigma # ======================================== # # Solve Eq. (17) of Li & Ma (1983) for Non # # ======================================== # def Non_lima(sigma, Noff, alpha=0.2): - """ - Solve Eq. (17) of Li & Ma (1983) for Non - - Parameters: - sigma - Required significance - Noff - Number of off counts - Keywords: - alpha - Ratio of on-to-off exposure - """ - # Set initial guess for Non - S = 25.0 - Non = S + alpha*Noff - - # Compute initial sigma estimate - s = sigma_lima(Non, Noff, alpha=alpha) - delta = s - sigma - - # Iterate - while abs(delta) > 1.0e-6: - S *= sigma/s - Non = S + alpha*Noff - s = sigma_lima(Non, Noff, alpha=alpha) - delta = s - sigma - - # Return - return Non + """ + Solve Eq. (17) of Li & Ma (1983) for Non + + Parameters: + sigma - Required significance + Noff - Number of off counts + Keywords: + alpha - Ratio of on-to-off exposure + """ + # Set initial guess for Non + S = 25.0 + Non = S + alpha*Noff + + # Compute initial sigma estimate + s = sigma_lima(Non, Noff, alpha=alpha) + delta = s - sigma + + # Iterate + while abs(delta) > 1.0e-6: + S *= sigma/s + Non = S + alpha*Noff + s = sigma_lima(Non, Noff, alpha=alpha) + delta = s - sigma + + # Return + return Non # =============================================== # @@ -97,195 +97,195 @@ def Non_lima(sigma, Noff, alpha=0.2): # =============================================== # def get_aeff(bgrate, emin, emax, \ duration=180000.0, alpha=0.2, sigma=5.0, norm=2.79e-7, index=-2.57): - """ - Computes effective area assuming a Crab spectrum. This does not - include any hadron cut efficiencies, and thus is not very useful in a first - place ... - - Parameters: - bgrate - Background rate (TeV) - emin - Minimum energy (TeV) - emax - Maximum energy (TeV) - Keywords: - duration - Observing time in seconds (default: 50h) - alpha - Non/Noff ratio (default: 0.2) - sigma - Significance (default: 5.0) - norm - Crab spectrum normalization at 1 TeV (default: 2.79e-7) - index - Crab spectral index (default: -2.57) - """ - # Compute number of source events - B = bgrate*duration - Noff = B/alpha - Non = Non_lima(sigma, Noff) - S = (Non - B) - - # Compute Crab flux in specified energy band - inx = index+1.0 - cu = norm/inx * (math.pow(emax, inx) - math.pow(emin, inx)) - - # Compute effective area - area = S / duration / cu - - # Return effective area - return area + """ + Computes effective area assuming a Crab spectrum. This does not + include any hadron cut efficiencies, and thus is not very useful in a first + place ... + + Parameters: + bgrate - Background rate (TeV) + emin - Minimum energy (TeV) + emax - Maximum energy (TeV) + Keywords: + duration - Observing time in seconds (default: 50h) + alpha - Non/Noff ratio (default: 0.2) + sigma - Significance (default: 5.0) + norm - Crab spectrum normalization at 1 TeV (default: 2.79e-7) + index - Crab spectral index (default: -2.57) + """ + # Compute number of source events + B = bgrate*duration + Noff = B/alpha + Non = Non_lima(sigma, Noff) + S = (Non - B) + + # Compute Crab flux in specified energy band + inx = index+1.0 + cu = norm/inx * (math.pow(emax, inx) - math.pow(emin, inx)) + + # Compute effective area + area = S / duration / cu + + # Return effective area + return area # ============== # # Generate files # # ============== # def make_files(filename, prfname, bkgname): - """ - Generate files. - - Parameters: - filename - Performance file name - prfname - Performance file name - bkgname - File function filename - """ - # Open performance file and file function - file = TFile(filename) - fprf = open(prfname, "w") - fbkg = open(bkgname, "w") - - # Write header - fprf.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n") - #sys.stdout.write("log(E) Area Area80 r68 r80 ERes. BG Rate Diff Sens BGRateSqDeg BGInt (BGIntControl)\n") + """ + Generate files. + + Parameters: + filename - Performance file name + prfname - Performance file name + bkgname - File function filename + """ + # Open performance file and file function + file = TFile(filename) + fprf = open(prfname, "w") + fbkg = open(bkgname, "w") + + # Write header + fprf.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n") + #sys.stdout.write("log(E) Area Area80 r68 r80 ERes. BG Rate Diff Sens BGRateSqDeg BGInt (BGIntControl)\n") - # Get histograms. Konrad's root files do not have an EffectiveArea80 - # histogram, but his effective area is for 80% containment radius - # anyways, so we simply read this histogram into aeff80. - sens = TH1F() - bgrate = TH1F() - bgratesqdeg = TH1F() - aeff = TH1F() - aeff80 = TH1F() - angres = TH1F() - angres80 = TH1F() - eres = TH1F() - file.GetObject("DiffSens", sens) - file.GetObject("BGRate", bgrate) - file.GetObject("BGRatePerSqDeg", bgratesqdeg) - file.GetObject("EffectiveArea", aeff) - try: - file.GetObject("EffectiveArea80", aeff80) - except: - file.GetObject("EffectiveArea", aeff80) - file.GetObject("AngRes", angres) - file.GetObject("AngRes80", angres80) - file.GetObject("ERes", eres) - - # Setup vectors - energies = bgrate.GetXaxis() - nbins_eng = energies.GetNbins() - for i in range(nbins_eng): - - # Get logE and background rate per squared degree - logE = energies.GetBinCenter(i+1) - bgd = bgrate.GetBinContent(i+1,1) - bgdsq = bgratesqdeg.GetBinContent(i+1,1) - area = aeff.GetBinContent(i+1,1) - area80 = aeff80.GetBinContent(i+1,1) - r68 = angres.GetBinContent(i+1,1) - r80 = angres80.GetBinContent(i+1,1) - dE = eres.GetBinContent(i+1,1) - diffs = sens.GetBinContent(i+1,1) - - # Skip any NaNs - if math.isnan(area80): - continue + # Get histograms. Konrad's root files do not have an EffectiveArea80 + # histogram, but his effective area is for 80% containment radius + # anyways, so we simply read this histogram into aeff80. + sens = TH1F() + bgrate = TH1F() + bgratesqdeg = TH1F() + aeff = TH1F() + aeff80 = TH1F() + angres = TH1F() + angres80 = TH1F() + eres = TH1F() + file.GetObject("DiffSens", sens) + file.GetObject("BGRate", bgrate) + file.GetObject("BGRatePerSqDeg", bgratesqdeg) + file.GetObject("EffectiveArea", aeff) + try: + file.GetObject("EffectiveArea80", aeff80) + except: + file.GetObject("EffectiveArea", aeff80) + file.GetObject("AngRes", angres) + file.GetObject("AngRes80", angres80) + file.GetObject("ERes", eres) + + # Setup vectors + energies = bgrate.GetXaxis() + nbins_eng = energies.GetNbins() + for i in range(nbins_eng): + + # Get logE and background rate per squared degree + logE = energies.GetBinCenter(i+1) + bgd = bgrate.GetBinContent(i+1,1) + bgdsq = bgratesqdeg.GetBinContent(i+1,1) + area = aeff.GetBinContent(i+1,1) + area80 = aeff80.GetBinContent(i+1,1) + r68 = angres.GetBinContent(i+1,1) + r80 = angres80.GetBinContent(i+1,1) + dE = eres.GetBinContent(i+1,1) + diffs = sens.GetBinContent(i+1,1) + + # Skip any NaNs + if math.isnan(area80): + continue - # Compute energy (in MeV) - energy = math.pow(10.0, logE)*1.0e6 - emin = math.pow(10.0, logE-0.1)*1.0e6 - emax = math.pow(10.0, logE+0.1)*1.0e6 - ewidth = emax-emin - - # Convert into background rate per steradian and MeV - bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth - - # Compute control background rate (this only works for KB files - # as they are for 80% containment radius). But we don't use this - # anyways, it's just for control on display ... - omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80))) - if omega > 0.0: - bkg_rate2 = bgd / omega / ewidth - else: - bkg_rate2 = 0.0 - - # Compute full effective area by dividing area80/0.80 - area = area80/0.80 + # Compute energy (in MeV) + energy = math.pow(10.0, logE)*1.0e6 + emin = math.pow(10.0, logE-0.1)*1.0e6 + emax = math.pow(10.0, logE+0.1)*1.0e6 + ewidth = emax-emin + + # Convert into background rate per steradian and MeV + bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth + + # Compute control background rate (this only works for KB files + # as they are for 80% containment radius). But we don't use this + # anyways, it's just for control on display ... + omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80))) + if omega > 0.0: + bkg_rate2 = bgd / omega / ewidth + else: + bkg_rate2 = 0.0 + + # Compute full effective area by dividing area80/0.80 + area = area80/0.80 - # Write results in file - line = "%.1f %.1f %.4f %.4f %.4f %.5e %.5e" % \ - (logE, area, r68, r80, dE, bgd, diffs) - fprf.write(line+"\n") - fbkg.write(str(energy)+" "+str(bkg_rate)+"\n") + # Write results in file + line = "%.1f %.1f %.4f %.4f %.4f %.5e %.5e" % \ + (logE, area, r68, r80, dE, bgd, diffs) + fprf.write(line+"\n") + fbkg.write(str(energy)+" "+str(bkg_rate)+"\n") - # Show performance file - #print "%.1f %.1f %.1f %.4f %.4f %.4f %.7f %.5e %.5e %.5e (%.5e)" % \ - # (logE, area, area80, r68, r80, dE, bgd, diffs, bgdsq, bkg_rate, bkg_rate2) + # Show performance file + #print "%.1f %.1f %.1f %.4f %.4f %.4f %.7f %.5e %.5e %.5e (%.5e)" % \ + # (logE, area, area80, r68, r80, dE, bgd, diffs, bgdsq, bkg_rate, bkg_rate2) - # Write trailer - fprf.write("---------------------------------------------\n") - fprf.write("Notes\n") - fprf.write(" 1) log(E) = log10(E/TeV) - bin centre\n") - fprf.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n") - fprf.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n") - fprf.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n") - fprf.write(" 5) Fractional Energy Resolution (rms)\n") - fprf.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n") - fprf.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n") - fprf.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n") - fprf.write(" systematics and statistics and at least 10 photons.\n") + # Write trailer + fprf.write("---------------------------------------------\n") + fprf.write("Notes\n") + fprf.write(" 1) log(E) = log10(E/TeV) - bin centre\n") + fprf.write(" 2) Eff Area - in square metres after background cut (no theta cut)\n") + fprf.write(" 3) Ang. Res - 68% containment radius of gamma-ray PSF post cuts - in degrees\n") + fprf.write(" 4) Ang. Res - 80% containment radius of gamma-ray PSF post cuts - in degrees\n") + fprf.write(" 5) Fractional Energy Resolution (rms)\n") + fprf.write(" 6) BG Rate - inside point-source selection region - post call cuts - in Hz\n") + fprf.write(" 7) Diff Sens - differential sensitivity for this bin expressed as E^2 dN/dE\n") + fprf.write(" - in erg cm^-2 s^-1 - for a 50 hours exposure - 5 sigma significance including\n") + fprf.write(" systematics and statistics and at least 10 photons.\n") - # Close files - fprf.close() - fbkg.close() + # Close files + fprf.close() + fbkg.close() - # Return - return + # Return + return #==========================# # Main routine entry point # #==========================# if __name__ == '__main__': - """ - """ - # Get ROOT filename from command line argument - if (len(sys.argv) > 1): - filename = sys.argv[1] - else: - #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEPerformanceBCE_Nov2010" - #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_subarrays_root" - caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_50deg_root" + """ + """ + # Get ROOT filename from command line argument + if (len(sys.argv) > 1): + filename = sys.argv[1] + else: + #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/IFAEPerformanceBCE_Nov2010" + #caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_subarrays_root" + caldb = "/Users/jurgen/Documents/Travail/Projects/CTA/WP-MC/root/Kb_all_50deg_root" - # Get list of all performance files - files = glob.glob(caldb+"/*.root") - files.sort() - - # Loop over all files - for file in files: - - # Get filename - head, tail = os.path.split(file) - - # Build name - name = tail.strip(".root") - if name.find("kb_") != -1: - name = name.replace("_20deg", "") - elif name.find("IFAE_") != -1: - name = name.replace("Subarray", "") - name = name.replace("IFAE_", "") - name = "IFAE_"+name - i = name.find("hours") - if i != -1: - name = name[0:i+1] - - # Build output filenames - prfname = name + ".dat" - bkgname = "bkg_"+name+".txt" - #print prfname, bkgname + # Get list of all performance files + files = glob.glob(caldb+"/*.root") + files.sort() + + # Loop over all files + for file in files: + + # Get filename + head, tail = os.path.split(file) + + # Build name + name = tail.strip(".root") + if name.find("kb_") != -1: + name = name.replace("_20deg", "") + elif name.find("IFAE_") != -1: + name = name.replace("Subarray", "") + name = name.replace("IFAE_", "") + name = "IFAE_"+name + i = name.find("hours") + if i != -1: + name = name[0:i+1] + + # Build output filenames + prfname = name + ".dat" + bkgname = "bkg_"+name+".txt" + #print prfname, bkgname - # Create files - make_files(file, prfname, bkgname) + # Create files + make_files(file, prfname, bkgname) diff --git a/src/ctlike/ctlike.cpp b/src/ctlike/ctlike.cpp index d003648..f9e3dff 100644 --- a/src/ctlike/ctlike.cpp +++ b/src/ctlike/ctlike.cpp @@ -356,15 +356,18 @@ void ctlike::get_parameters(void) // Get event file name std::string filename = (*this)["infile"].filename(); - // Try first to open as unbinned FITS file + // Try first to open as FITS file try { // Determine whether FITS file has EVENTS extension + /* GFits fits(filename); bool is_unbinned = fits.contains("EVENTS"); fits.close(); + */ // If FITS file has an events header then load as unbinned + /* if (is_unbinned) { obs.load_unbinned(filename); } @@ -373,8 +376,12 @@ void ctlike::get_parameters(void) else { obs.load_binned(filename); } + */ - // Get task parameters for event file loading + // Load data + obs.load(filename); + + // Get other task parameters m_stat = gammalib::toupper((*this)["stat"].string()); m_caldb = (*this)["caldb"].string(); m_irf = (*this)["irf"].string(); @@ -408,7 +415,45 @@ void ctlike::get_parameters(void) // Load observations from XML file m_obs.load(filename); - } + // Check if all observations have response information. If + // not, get the calibration database parameters and set + // the response properly + bool asked_for_response = false; + for (int i = 0; i < m_obs.size(); ++i) { + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(m_obs[i]); + + // Continue only if observation is a CTA observation + if (obs != NULL) { + + // If there is no effective area information then + // set response from task parameters + if (obs->response().aeff() == NULL) { + + // Get calibration parameters + if (!asked_for_response) { + m_caldb = (*this)["caldb"].string(); + m_irf = (*this)["irf"].string(); + asked_for_response = true; + } + + // Set response + GCaldb caldb; + if (gammalib::dir_exists(m_caldb)) { + caldb.rootdir(m_caldb); + } + else { + caldb.open("cta", m_caldb); + } + obs->response(m_irf, caldb); + + } // endif: there was no response + } // endif: observation was a CTA observation + + } // endfor: looped over observations + + } // endcatch: file was an XML file } // endif: there was no observation in the container diff --git a/src/ctlike/ctlike.par b/src/ctlike/ctlike.par index 85d939b..cdaf144 100644 --- a/src/ctlike/ctlike.par +++ b/src/ctlike/ctlike.par @@ -17,14 +17,13 @@ # # General parameters #=================== -infile, f, a, "events.fits",,,"Event list, counts map or observation definition file" +infile, f, a, "events.fits",,, "Event list, counts map or observation definition file" stat, s, h, "POISSON",,, "Optimization statistics" refit, b, h, no,,, "Do refitting?" -edisp, b, h, no,,, "Apply energy dispersion?" -caldb, s, a, "$CTOOLS/share/caldb/data/cta/dummy",,, "Calibration database" +caldb, s, a, "dummy",,, "Calibration database" irf, s, a, "cta_dummy_irf",,, "Instrument response function" srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" -outmdl, f, a, "crab_results.xml",,, "Source model output file" +outmdl, f, a, "crab_results.xml",,, "Source model output file" # # Standard parameters #==================== diff --git a/src/ctobssim/ctobssim.par b/src/ctobssim/ctobssim.par index c246ccf..19e9e49 100644 --- a/src/ctobssim/ctobssim.par +++ b/src/ctobssim/ctobssim.par @@ -18,12 +18,12 @@ # File information #================= infile, f, a, "$CTOOLS/share/models/crab.xml",,, "Model" -outfile,f, a, "events.fits",,, "Output event data file or observation definition file" -prefix, s, h, "sim_events_",,, "Prefix for event lists in observation definition file" +outfile,f, a, "events.fits",,, "Output event data file or observation definition file" +prefix, s, h, "sim_events_",,, "Prefix for event lists in observation definition file" # # Observation parameters #======================= -caldb, s, a, "$CTOOLS/share/caldb/data/cta/dummy",,, "Calibration database" +caldb, s, a, "dummy",,, "Calibration database" irf, s, a, "cta_dummy_irf",,, "Instrument response function" # # Simulation parameters @@ -31,7 +31,7 @@ irf, s, a, "cta_dummy_irf",,, "Instrument response function" seed, i, h, 1, , , "Random number generator seed" ra, r, a, 83.63, 0,360, "RA of pointing (degrees)" dec, r, a, 22.01, -90, 90, "Dec of pointing (degrees)" -rad, r, a, 10.0, 0,180, "Radius of FOV (degrees)" +rad, r, a, 5.0, 0,180, "Radius of FOV (degrees)" tmin, r, a, 0.0, 0, , "Start time (MET in s)" tmax, r, a, 1800.0, 0, , "End time (MET in s)" emin, r, a, 0.1, 0, , "Lower energy limit (TeV)"