diff --git a/doc/source/users/tutorials/hess_dr1/fitting.rst b/doc/source/users/tutorials/hess_dr1/fitting.rst new file mode 100644 index 0000000..6866633 --- /dev/null +++ b/doc/source/users/tutorials/hess_dr1/fitting.rst @@ -0,0 +1,95 @@ +.. _hess_dr1_fitting: + +Fitting the model components to the data +---------------------------------------- + + .. admonition:: What you will learn + + You will learn **how to fit a parametric model to the data**. + +Adjust source model name and set pivot energy to 1 TeV to bring it energy +interval of data. +Remove ``error`` attributes. + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + + + + + +.. code-block:: bash + + $ ctlike debug=yes + Input event list, counts cube or observation definition XML file [events.fits] obs_crab_selected.xml + Input model definition XML file [$CTOOLS/share/models/crab.xml] crab_models.xml + Output model definition XML file [crab_results.xml] + +.. code-block:: none + + 2018-09-11T21:17:31: +=========================================+ + 2018-09-11T21:17:31: | Maximum likelihood optimisation results | + 2018-09-11T21:17:31: +=========================================+ + 2018-09-11T21:17:31: === GOptimizerLM === + 2018-09-11T21:17:31: Optimized function value ..: 98422.688 + 2018-09-11T21:17:31: Absolute precision ........: 0.005 + 2018-09-11T21:17:31: Acceptable value decrease .: 2 + 2018-09-11T21:17:31: Optimization status .......: converged + 2018-09-11T21:17:31: Number of parameters ......: 10 + 2018-09-11T21:17:31: Number of free parameters .: 6 + 2018-09-11T21:17:31: Number of iterations ......: 24 + 2018-09-11T21:17:31: Lambda ....................: 0.1 + 2018-09-11T21:17:31: Maximum log likelihood ....: -98422.688 + 2018-09-11T21:17:31: Observed events (Nobs) ...: 9675.000 + 2018-09-11T21:17:31: Predicted events (Npred) ..: 9674.870 (Nobs - Npred = 0.129601074213497) + 2018-09-11T21:17:31: === GModels === + 2018-09-11T21:17:31: Number of models ..........: 2 + 2018-09-11T21:17:31: Number of parameters ......: 10 + 2018-09-11T21:17:31: === GModelSky === + 2018-09-11T21:17:31: Name ......................: Crab + 2018-09-11T21:17:31: Instruments ...............: all + 2018-09-11T21:17:31: Instrument scale factors ..: unity + 2018-09-11T21:17:31: Observation identifiers ...: all + 2018-09-11T21:17:31: Model type ................: PointSource + 2018-09-11T21:17:31: Model components ..........: "PointSource" * "PowerLaw" * "Constant" + 2018-09-11T21:17:31: Number of parameters ......: 6 + 2018-09-11T21:17:31: Number of spatial par's ...: 2 + 2018-09-11T21:17:31: RA .......................: 83.6230374464172 +/- 0.00241074887155303 deg (free,scale=1) + 2018-09-11T21:17:31: DEC ......................: 22.0249049416179 +/- 0.0022137104783551 deg (free,scale=1) + 2018-09-11T21:17:31: Number of spectral par's ..: 3 + 2018-09-11T21:17:31: Prefactor ................: 4.81375364083447e-17 +/- 2.62970930530696e-18 [0,infty[ ph/cm2/s/MeV (free,scale=5.7e-18,gradient) + 2018-09-11T21:17:31: Index ....................: -2.71079050853091 +/- 0.0652734279744327 [9.92,-9.92] (free,scale=-2.48,gradient) + 2018-09-11T21:17:31: PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient) + 2018-09-11T21:17:31: Number of temporal par's ..: 1 + 2018-09-11T21:17:31: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2018-09-11T21:17:31: === GCTAModelAeffBackground === + 2018-09-11T21:17:31: Name ......................: Background + 2018-09-11T21:17:31: Instruments ...............: all + 2018-09-11T21:17:31: Instrument scale factors ..: unity + 2018-09-11T21:17:31: Observation identifiers ...: all + 2018-09-11T21:17:31: Model type ................: "PowerLaw" * "Constant" + 2018-09-11T21:17:31: Number of parameters ......: 4 + 2018-09-11T21:17:31: Number of spectral par's ..: 3 + 2018-09-11T21:17:31: Prefactor ................: 1.80502240642969e-13 +/- 2.34991048262393e-15 [0,infty[ ph/cm2/s/MeV (free,scale=1e-13,gradient) + 2018-09-11T21:17:31: Index ....................: -2.58303579414564 +/- 0.0128891124984818 [10,-10] (free,scale=-2.5,gradient) + 2018-09-11T21:17:31: PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient) + 2018-09-11T21:17:31: Number of temporal par's ..: 1 + 2018-09-11T21:17:31: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) diff --git a/doc/source/users/tutorials/hess_dr1/index.rst b/doc/source/users/tutorials/hess_dr1/index.rst index c29d0bf..740515e 100644 --- a/doc/source/users/tutorials/hess_dr1/index.rst +++ b/doc/source/users/tutorials/hess_dr1/index.rst @@ -7,4 +7,6 @@ Analysing public H.E.S.S. data :maxdepth: 1 data - + select_events + skymap + fitting diff --git a/doc/source/users/tutorials/hess_dr1/select_events.rst b/doc/source/users/tutorials/hess_dr1/select_events.rst new file mode 100644 index 0000000..7835a7b --- /dev/null +++ b/doc/source/users/tutorials/hess_dr1/select_events.rst @@ -0,0 +1,145 @@ +.. _hess_dr1_select_events: + +Selecting the relevant events +----------------------------- + + .. admonition:: What you will learn + + You will learn **how to select the relevant events for a target** + and **how to determine the useful energy range** for a set of + observations. + +Now you are ready for your first H.E.S.S. data analysis. +As first step you should create a folder that will contain your analysis +results and any intermediate data products. +Let's call this folder ``crab`` since we explain in this example how you +can analyse the data from the Crab nebula. +Create and step into the folder using + +.. code-block:: bash + + $ mkdir crab + $ cd crab + +The first step of your analysis consists in selecting the relevant events +from the observations. +In this step you can select a specific energy range, time range or region +of interest. +In the example below you will select events according to the save energy +thresholds that are defined in the effective area component of the +:ref:`instrument response functions `. +You do this by setting the hidden parameter ``usethres=DEFAULT``. +In addition, you will select all events within 2 degrees of the pointing +direction: + +.. code-block:: bash + + $ ctselect usethres=DEFAULT + Input event list or observation definition XML file [events.fits] $HESSDATA/obs/obs_crab.xml + Lower energy limit (TeV) [0.1] NONE + Radius of ROI around pointing or specified RA/DEC (degrees) (0-180) [3.0] 2.0 + Start time (UTC string, JD, MJD or MET in seconds) [NONE] + Output event list or observation definition XML file [selected_events.fits] obs_crab_selected.xml + +Below is an excerpt of the log file that was created by :ref:`ctselect`. +You will find that the following energy ranges were applied for the +observations: + +* ``id=023523``: 0.871 - 100 TeV +* ``id=023526``: 0.708 - 100 TeV +* ``id=023559``: 0.661 - 100 TeV +* ``id=023592``: 0.871 - 100 TeV + +.. code-block:: none + + 2018-09-11T14:05:15: +=================+ + 2018-09-11T14:05:15: | Event selection | + 2018-09-11T14:05:15: +=================+ + 2018-09-11T14:05:15: === HESS observation "Crab" (id=023523) === + 2018-09-11T14:05:15: Input filename ............: /project-data/hess/hess_dl3_dr1/data/hess_dl3_dr1_obs_id_023523.fits.gz + 2018-09-11T14:05:15: Event extension name ......: EVENTS + 2018-09-11T14:05:15: GTI extension name ........: GTI + 2018-09-11T14:05:16: === Event selection === + 2018-09-11T14:05:16: Selected energy range .....: 0.870964 - 100 TeV + 2018-09-11T14:05:16: Requested RoI .............: Centre(RA,DEC)=(83.633333333333, 21.514444444444) deg, Radius=2 deg + 2018-09-11T14:05:16: RoI of data ...............: Centre(RA,DEC)=(83.633333333333, 21.514444444444) deg, Radius=2 deg + 2018-09-11T14:05:16: Selected RoI ..............: Centre(RA,DEC)=(83.633333333333, 21.514444444444) deg, Radius=2 deg + 2018-09-11T14:05:16: cfitsio selection .........: ENERGY >= 0.87096359 && ENERGY <= 100.00000000 && ANGSEP(83.633333,21.514444,RA,DEC) <= 2.000000 + 2018-09-11T14:05:16: FITS filename .............: /var/tmp/tmp.0.QX6f0L[EVENTS][ENERGY >= 0.87096359 && ENERGY <= 100.00000000 && ANGSEP(83.633333,21.514444,RA,DEC) <= 2.000000] + 2018-09-11T14:05:16: === HESS observation "Crab" (id=023526) === + 2018-09-11T14:05:16: Input filename ............: /project-data/hess/hess_dl3_dr1/data/hess_dl3_dr1_obs_id_023526.fits.gz + 2018-09-11T14:05:16: Event extension name ......: EVENTS + 2018-09-11T14:05:16: GTI extension name ........: GTI + 2018-09-11T14:05:16: === Event selection === + 2018-09-11T14:05:16: Selected energy range .....: 0.707946 - 100 TeV + 2018-09-11T14:05:16: Requested RoI .............: Centre(RA,DEC)=(83.633333333333, 22.514444444444) deg, Radius=2 deg + 2018-09-11T14:05:16: RoI of data ...............: Centre(RA,DEC)=(83.633333333333, 22.514444444444) deg, Radius=2 deg + 2018-09-11T14:05:16: Selected RoI ..............: Centre(RA,DEC)=(83.633333333333, 22.514444444444) deg, Radius=2 deg + 2018-09-11T14:05:16: cfitsio selection .........: ENERGY >= 0.70794578 && ENERGY <= 100.00000000 && ANGSEP(83.633333,22.514444,RA,DEC) <= 2.000000 + 2018-09-11T14:05:16: FITS filename .............: /var/tmp/tmp.1.IoZj2b[EVENTS][ENERGY >= 0.70794578 && ENERGY <= 100.00000000 && ANGSEP(83.633333,22.514444,RA,DEC) <= 2.000000] + 2018-09-11T14:05:16: === HESS observation "Crab" (id=023559) === + 2018-09-11T14:05:16: Input filename ............: /project-data/hess/hess_dl3_dr1/data/hess_dl3_dr1_obs_id_023559.fits.gz + 2018-09-11T14:05:16: Event extension name ......: EVENTS + 2018-09-11T14:05:16: GTI extension name ........: GTI + 2018-09-11T14:05:16: === Event selection === + 2018-09-11T14:05:16: Selected energy range .....: 0.660693 - 100 TeV + 2018-09-11T14:05:16: Requested RoI .............: Centre(RA,DEC)=(85.2533333381014, 22.014444444444) deg, Radius=2 deg + 2018-09-11T14:05:16: RoI of data ...............: Centre(RA,DEC)=(85.2533333381014, 22.014444444444) deg, Radius=2 deg + 2018-09-11T14:05:16: Selected RoI ..............: Centre(RA,DEC)=(85.2533333381014, 22.014444444444) deg, Radius=2 deg + 2018-09-11T14:05:16: cfitsio selection .........: ENERGY >= 0.66069345 && ENERGY <= 100.00000000 && ANGSEP(85.253333,22.014444,RA,DEC) <= 2.000000 + 2018-09-11T14:05:16: FITS filename .............: /var/tmp/tmp.2.IAwBNH[EVENTS][ENERGY >= 0.66069345 && ENERGY <= 100.00000000 && ANGSEP(85.253333,22.014444,RA,DEC) <= 2.000000] + 2018-09-11T14:05:16: === HESS observation "Crab" (id=023592) === + 2018-09-11T14:05:16: Input filename ............: /project-data/hess/hess_dl3_dr1/data/hess_dl3_dr1_obs_id_023592.fits.gz + 2018-09-11T14:05:16: Event extension name ......: EVENTS + 2018-09-11T14:05:16: GTI extension name ........: GTI + 2018-09-11T14:05:16: === Event selection === + 2018-09-11T14:05:16: Selected energy range .....: 0.870964 - 100 TeV + 2018-09-11T14:05:16: Requested RoI .............: Centre(RA,DEC)=(82.0133333285646, 22.014444444444) deg, Radius=2 deg + 2018-09-11T14:05:16: RoI of data ...............: Centre(RA,DEC)=(82.0133333285646, 22.014444444444) deg, Radius=2 deg + 2018-09-11T14:05:16: Selected RoI ..............: Centre(RA,DEC)=(82.0133333285646, 22.014444444444) deg, Radius=2 deg + 2018-09-11T14:05:16: cfitsio selection .........: ENERGY >= 0.87096359 && ENERGY <= 100.00000000 && ANGSEP(82.013333,22.014444,RA,DEC) <= 2.000000 + 2018-09-11T14:05:16: FITS filename .............: /var/tmp/tmp.3.879wq4[EVENTS][ENERGY >= 0.87096359 && ENERGY <= 100.00000000 && ANGSEP(82.013333,22.014444,RA,DEC) <= 2.000000] + +You may use the :ref:`csobsinfo` script to display a summary of the selected +events into the console: + +.. code-block:: bash + + $ csobsinfo debug=yes + Input event list, counts cube, or observation definition XML file [$HESSDATA/obs/obs_crab.xml] obs_crab_selected.xml + Output DS9 region file [ds9.reg] + +.. note:: + The ``debug=yes`` attribute instructs :ref:`csobsinfo` to direct the log + file output also to the console. Duplication of log file output into + the console using the ``debug=yes`` attribute works for any tool or script. + +The output of :ref:`csobsinfo` is shown below. +:ref:`ctselect` selected 9675 events within an energy interval of +0.66 - 100 TeV. +This will be the energy range that you will from now on use for the data +analysis. + +.. code-block:: none + + 2018-09-11T20:30:16: +=========+ + 2018-09-11T20:30:16: | Summary | + 2018-09-11T20:30:16: +=========+ + 2018-09-11T20:30:16: === Observations === + 2018-09-11T20:30:16: Unbinned observations .....: 4 + 2018-09-11T20:30:16: Binned observations .......: 0 + 2018-09-11T20:30:16: === Events === + 2018-09-11T20:30:16: Number of events ..........: 9675 + 2018-09-11T20:30:16: Number of bins ............: 0 + 2018-09-11T20:30:16: === Pointings === + 2018-09-11T20:30:16: Mean offset angle .........: Unknown + 2018-09-11T20:30:16: Mean zenith angle .........: 47.04 deg + 2018-09-11T20:30:16: Mean azimuth angle ........: 13.76 deg + 2018-09-11T20:30:16: === Energy range === + 2018-09-11T20:30:16: Minimum energy ............: 660.693448007596 GeV + 2018-09-11T20:30:16: Maximum energy ............: 100 TeV + 2018-09-11T20:30:16: === Time range === + 2018-09-11T20:30:16: MJD (days) ................: 53343.922 - 53347.933 + 2018-09-11T20:30:16: UTC .......................: 2004-12-04T22:07:06 - 2004-12-08T22:22:02 + 2018-09-11T20:30:16: Total ontime ..............: 6742.00 s = 112.37 min = 1.87 h + 2018-09-11T20:30:16: Total livetime ............: 6313.81 s = 105.23 min = 1.75 h diff --git a/doc/source/users/tutorials/hess_dr1/skymap.png b/doc/source/users/tutorials/hess_dr1/skymap.png new file mode 100644 index 0000000..5b4bebc --- /dev/null +++ b/doc/source/users/tutorials/hess_dr1/skymap.png diff --git a/doc/source/users/tutorials/hess_dr1/skymap.rst b/doc/source/users/tutorials/hess_dr1/skymap.rst new file mode 100644 index 0000000..7b71ec4 --- /dev/null +++ b/doc/source/users/tutorials/hess_dr1/skymap.rst @@ -0,0 +1,124 @@ +.. _hess_dr1_sky_map: + +Generating a sky map from the events +------------------------------------ + + .. admonition:: What you will learn + + You will learn **how to generate a sky map from the selected events** + and **how to find potential sources in the data**. + +As next step you will make a sky map of the selected events to inspect what's +in the data. +You do this using the :ref:`ctskymap` tool. +Since the data correspond to observations of the Crab nebula you should centre +the sky map on this source. +The energy limits of the sky map should be set to the energy interval of the +data (0.66 - 100 TeV). +To create the sky map, type: + +.. code-block:: bash + + $ ctskymap debug=yes + Input event list or observation definition XML file [events.fits] obs_crab_selected.xml + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] + Projection method (AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN) [CAR] + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] + Image scale (in degrees/pixel) [0.02] + Size of the X axis in pixels [200] + Size of the Y axis in pixels [200] + Lower energy limit (TeV) [0.1] 0.66 + Upper energy limit (TeV) [100.0] + Background subtraction method (NONE|IRF|RING) [NONE] + Output skymap file [skymap.fits] + +The resulting sky map, displayed using +`ds9 `_ with a square-root color scaling, is shown below: + +.. figure:: skymap.png + :width: 400px + :align: center + + *Sky map of the selected events centred on the Crab nebula* + +Based on the sky map you can now extract a model for all potentially significant +sources in the data. +You do this using the :ref:`cssrcdetect` script, and to extract only significant +sources, you should set the detection threshold for example to 5 sigma. +For convenience, the :ref:`cssrcdetect` script can also add a background model +component to the output +:ref:`model definition XML file `. +Since the H.E.S.S. Instrument Response Function does not include a background +template, the ``IRF`` option will not work. +While sophisticated methods for the analysis of H.E.S.S. data still need to +be developed, you may try the effective area model ``AEFF`` that uses the +effective area as spatial template for the background component. +So, now type: + +.. code-block:: bash + + $ cssrcdetect + Input sky map file [skymap.fits] + Source model type (POINT) [POINT] + Background model type (NONE|IRF|AEFF|CUBE|RACC) [NONE] AEFF + Detection threshold (Gaussian sigma) [5.0] + Correlation kernel radius (deg) [0.1] + Output model definition XML file [models.xml] crab_models.xml + Output DS9 region file [ds9.reg] crab.reg + +Below an excerpt of the ``cssrcdetect.log`` log file. +Only one significant source was found, and its position is consistent with the +Crab nebula. +This is also illustrated in the figure below where the region file +``crab.reg`` that was produced by :ref:`cssrcdetect` was overlaid with +`ds9 `_ +on the sky map. + +.. code-block:: none + + 2018-09-11T20:32:20: +==================+ + 2018-09-11T20:32:20: | Source detection | + 2018-09-11T20:32:20: +==================+ + 2018-09-11T20:32:20: === Iteration 1 === + 2018-09-11T20:32:20: Map maximum ...............: 11.1003120131 + 2018-09-11T20:32:20: Src001 position ...........: (RA,Dec)=(83.6192131308071,22.0199996472185) + 2018-09-11T20:32:20: === Iteration 2 === + 2018-09-11T20:32:20: Map maximum ...............: None above threshold + +.. figure:: skymap_crab.png + :width: 400px + :align: center + + *Sky map of the selected events centred on the Crab nebula* + +The +:ref:`model definition file ` +``crab_models.xml`` that was generated by :ref:`cssrcdetect` is shown below. +It contains a point source ``Src001`` at the position of the Crab nebula +and an effective area background model scaled by a spectral power law as +background component. + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/source/users/tutorials/hess_dr1/skymap_crab.png b/doc/source/users/tutorials/hess_dr1/skymap_crab.png new file mode 100644 index 0000000..08c3fb4 --- /dev/null +++ b/doc/source/users/tutorials/hess_dr1/skymap_crab.png