diff --git a/cscripts/csiactobs.par b/cscripts/csiactobs.par index 5c3f7bd..073d18c 100755 --- a/cscripts/csiactobs.par +++ b/cscripts/csiactobs.par @@ -22,7 +22,7 @@ infile, f, a, runlist.lis,,, "Input runlist file" inmodel, f, h, NONE,,, "Input model definition XML file (optional)" master_indx, s, h, master.json,,, "Input master index file" outobs, f, a, obs.xml,,, "Output observation definition XML file" -outmodel, f, a, bgmodels.xml,,, "Output model definition XML file" +outmodel, f, a, bkgmodels.xml,,, "Output model definition XML file" # # Script parameters @@ -31,7 +31,7 @@ bkgpars, i, a, 1,,, "Number of free parameters per background model" bkg_scale, b, h, yes,,, "Apply background scaling factor from observation index file if available?" ev_hiera, s, h, events,,, "Hierarchy of event formats" aeff_hiera, s, h, aeff_2d,,, "Hierarchy of effective area formats" -psf_hiera, s, h, psf_king|psf_3gauss,,, "Hierarchy of psf formats" +psf_hiera, s, h, psf_king|psf_3gauss|psf_table,,, "Hierarchy of psf formats" edisp_hiera, s, h, edisp_2d,,, "Hierarchy of energy dispersion formats" bkg_hiera, s, h, bkg_3d,,, "Hierarchy of background formats" bkg_mod_hiera, s, h, irf|aeff|gauss,,, "Hierarchy of background models" diff --git a/doc/source/users/tutorials/hess_dr1/index.rst b/doc/source/users/tutorials/hess_dr1/index.rst index 6330107..e228b07 100644 --- a/doc/source/users/tutorials/hess_dr1/index.rst +++ b/doc/source/users/tutorials/hess_dr1/index.rst @@ -1,4 +1,4 @@ -.. _1dc: +.. _hess_dr1: Analysing public H.E.S.S. data ============================== diff --git a/doc/source/users/tutorials/hess_dr1/jointfermi.rst b/doc/source/users/tutorials/hess_dr1/jointfermi.rst index 511246c..63672ab 100644 --- a/doc/source/users/tutorials/hess_dr1/jointfermi.rst +++ b/doc/source/users/tutorials/hess_dr1/jointfermi.rst @@ -1,4 +1,4 @@ -.. _hess_dr1_data: +.. _hess_dr1_fermi: Joint analysis of H.E.S.S. and Fermi data ========================================= diff --git a/doc/source/users/tutorials/howto/iact/analysis.rst b/doc/source/users/tutorials/howto/iact/analysis.rst index 8b571b2..dd50e1e 100755 --- a/doc/source/users/tutorials/howto/iact/analysis.rst +++ b/doc/source/users/tutorials/howto/iact/analysis.rst @@ -1,572 +1,164 @@ .. _sec_iact_analysis: -Analysing IACT data -=================== - -In order to analyse data from current IACTs, you first have to retrieve the -data. Information on how to do that can be found in the section :ref:`sec_iact_copy`. - -Check available FITS production -------------------------------- -Before you start an analysis, it is important that you know what kind of data -need to be analysed. In the current scheme of -`storing IACT data `_, -every FITS data set is represented by a unique name, e.g. -``"hess-fits-pa-release-1.0-Prod26-MppStd"``. -You have to pass this string to some set of scripts to find and use the data -for analysis. - -To check the available names of FITS data sets, you should run the -:ref:`csiactdata` script: - -.. code-block:: bash - - $ csiactdata - Path were data is located [] /path/to/fits/data/ - -This script will list the names of available FITS data sets into the console. -In the following it is assumed that you have selected a data set named -``fits-data-name``. -You can omit the query for the path where the data is located by setting the -environment variable ``VHEFITS`` to this location. Also other scripts will -benefit from having this environment variable set, and this avoids that you -have to type in the path name every time you use one of these scripts. - - -Find observations ------------------ - -Now that you know which data are available you have to assemble a list of -observations (a.k.a. runs) that you want to use in an analysis. You will -use the :ref:`csfindobs` for that task, which allows you to select the -observations according to a number of user criteria. - - -Select observations according to pointing direction -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -The most common search for observations is by pointing direction on the sky, -and you do this by simply invoking the :ref:`csfindobs` script: - -.. code-block:: bash - - $ csfindobs - Name of FITS production (Run csiactdata to view your options) [fits-data-name] - Right ascension of selection region centre (deg) [83.63] - Declination of selection region centre (deg) [22.01] - Search radius of selection region (deg) [2.5] - Output runlist file [runlist.lis] - -The script will create an ASCII file named ``runlist.lis`` that contains the -identifiers for all observations with pointing directions located within 2.5° -around the position of the Crab nebula. - - -Select observations using additional criteria -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -You may add further selection criteria using the hidden ``expression`` parameter. -You may use any FITS selection expression that is -`supported by the cfitsio row filtering specification `_. -Column names that can be used in the expression can be found -`here `_ -or by simply browsing the observation index file with a FITS viewer. For example - -.. code-block:: bash - - $ csfindobs expression="ZEN_PNT<30&&LIVETIME>1800" - -will select all observations that have a zenith angle less than 30° and a -livetime larger than 1800 seconds. - - -Select observations independent of pointing direction -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -You can also omit the pointing selection, for example in case that you do not -know the sky coordinates of your object of interest. In that case, you may -directly use the object name in the expression and pass ``NONE`` as coordinate -to :ref:`csfindobs`: - -.. code-block:: bash - - $ csfindobs expression="OBJECT=='Crab Nebula'" - Name of FITS production (Run csiactdata to view your options) [fits-data-name] - Right ascension of selection region centre [NONE] - Output runlist file [runlist.lis] +Performing a maximum likelihood fit +=================================== + + .. admonition:: What you will learn + + You will learn how to perform a maximum likelihood fit of IACT data. + +Following a :ref:`csiactobs` run with the ``inmodel=$CTOOLS/share/models/crab.xml`` +you have a model that is ready for maximum likelihood fitting. +You do such a fit using + +.. code-block:: bash + + $ ctlike + Input event list, counts cube or observation definition XML file [obs_selected.xml] + Input model definition XML file [models.xml] + Output model definition XML file [crab_results.xml] + +which will produce the following output in the log file ``ctlike.log``: + +.. code-block:: none + + 2019-04-05T14:18:17: +=================================+ + 2019-04-05T14:18:17: | Maximum likelihood optimisation | + 2019-04-05T14:18:17: +=================================+ + 2019-04-05T14:18:17: Parameter "Prefactor" has zero curvature. Fix parameter. + 2019-04-05T14:18:17: Parameter "Index" has zero curvature. Fix parameter. + 2019-04-05T14:18:17: >Iteration 0: -logL=203790.772, Lambda=1.0e-03 + 2019-04-05T14:18:17: >Iteration 1: -logL=192855.522, Lambda=1.0e-03, delta=10935.250, step=1.0e+00, max(|grad|)=-2040.622705 [Prefactor:0] + 2019-04-05T14:18:17: >Iteration 2: -logL=182999.110, Lambda=1.0e-04, delta=9856.412, step=1.0e+00, max(|grad|)=-981.683624 [Prefactor:0] + 2019-04-05T14:18:17: >Iteration 3: -logL=174691.964, Lambda=1.0e-05, delta=8307.146, step=1.0e+00, max(|grad|)=620.662474 [Index:9] + 2019-04-05T14:18:17: >Iteration 4: -logL=168496.493, Lambda=1.0e-06, delta=6195.471, step=1.0e+00, max(|grad|)=587.098865 [Index:9] + 2019-04-05T14:18:17: >Iteration 5: -logL=164898.487, Lambda=1.0e-07, delta=3598.006, step=1.0e+00, max(|grad|)=423.800306 [Index:5] + 2019-04-05T14:18:18: >Iteration 6: -logL=163621.925, Lambda=1.0e-08, delta=1276.562, step=1.0e+00, max(|grad|)=237.546032 [Index:1] + 2019-04-05T14:18:18: >Iteration 7: -logL=163423.407, Lambda=1.0e-09, delta=198.518, step=1.0e+00, max(|grad|)=102.646722 [Index:1] + 2019-04-05T14:18:18: >Iteration 8: -logL=163409.717, Lambda=1.0e-10, delta=13.690, step=1.0e+00, max(|grad|)=35.371246 [Index:1] + 2019-04-05T14:18:18: >Iteration 9: -logL=163408.542, Lambda=1.0e-11, delta=1.175, step=1.0e+00, max(|grad|)=11.708402 [Index:1] + 2019-04-05T14:18:18: >Iteration 10: -logL=163408.416, Lambda=1.0e-12, delta=0.126, step=1.0e+00, max(|grad|)=3.907383 [Index:1] + 2019-04-05T14:18:18: >Iteration 11: -logL=163408.402, Lambda=1.0e-13, delta=0.014, step=1.0e+00, max(|grad|)=1.309476 [Index:1] + 2019-04-05T14:18:18: >Iteration 12: -logL=163408.400, Lambda=1.0e-14, delta=0.002, step=1.0e+00, max(|grad|)=0.439462 [Index:1] + 2019-04-05T14:18:18: Free parameter "Prefactor" after convergence was reached with frozen parameter. + 2019-04-05T14:18:18: Free parameter "Index" after convergence was reached with frozen parameter. + 2019-04-05T14:18:18: + 2019-04-05T14:18:18: +=========================================+ + 2019-04-05T14:18:18: | Maximum likelihood optimisation results | + 2019-04-05T14:18:18: +=========================================+ + 2019-04-05T14:18:18: === GOptimizerLM === + 2019-04-05T14:18:18: Optimized function value ..: 163408.400 + 2019-04-05T14:18:18: Absolute precision ........: 0.005 + 2019-04-05T14:18:18: Acceptable value decrease .: 2 + 2019-04-05T14:18:18: Optimization status .......: converged + 2019-04-05T14:18:18: Number of parameters ......: 26 + 2019-04-05T14:18:18: Number of free parameters .: 12 + 2019-04-05T14:18:18: Number of iterations ......: 12 + 2019-04-05T14:18:18: Lambda ....................: 1e-15 + 2019-04-05T14:18:18: Maximum log likelihood ....: -163408.400 + 2019-04-05T14:18:18: Observed events (Nobs) ...: 18641.000 + 2019-04-05T14:18:18: Predicted events (Npred) ..: 18640.998 (Nobs - Npred = 0.00164345248049358) + 2019-04-05T14:18:18: === GModels === + 2019-04-05T14:18:18: Number of models ..........: 6 + 2019-04-05T14:18:18: Number of parameters ......: 26 + 2019-04-05T14:18:18: === GCTAModelAeffBackground === + 2019-04-05T14:18:18: Name ......................: bkg_23523 + 2019-04-05T14:18:18: Instruments ...............: HESS + 2019-04-05T14:18:18: Instrument scale factors ..: unity + 2019-04-05T14:18:18: Observation identifiers ...: 23523 + 2019-04-05T14:18:18: Model type ................: "PowerLaw" * "Constant" + 2019-04-05T14:18:18: Number of parameters ......: 4 + 2019-04-05T14:18:18: Number of spectral par's ..: 3 + 2019-04-05T14:18:18: Prefactor ................: 2.98790819794319e-13 +/- 4.43400659819813e-15 [1e-16,1e-12] ph/cm2/s/MeV (free,scale=1e-14,gradient) + 2019-04-05T14:18:18: Index ....................: -3.22702359177827 +/- 0.0221521174108544 [-5,5] (free,scale=1,gradient) + 2019-04-05T14:18:18: PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient) + 2019-04-05T14:18:18: Number of temporal par's ..: 1 + 2019-04-05T14:18:18: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2019-04-05T14:18:18: === GCTAModelAeffBackground === + 2019-04-05T14:18:18: Name ......................: bkg_23526 + 2019-04-05T14:18:18: Instruments ...............: HESS + 2019-04-05T14:18:18: Instrument scale factors ..: unity + 2019-04-05T14:18:18: Observation identifiers ...: 23526 + 2019-04-05T14:18:18: Model type ................: "PowerLaw" * "Constant" + 2019-04-05T14:18:18: Number of parameters ......: 4 + 2019-04-05T14:18:18: Number of spectral par's ..: 3 + 2019-04-05T14:18:18: Prefactor ................: 2.44289795381698e-13 +/- 3.75318255643637e-15 [1e-16,1e-12] ph/cm2/s/MeV (free,scale=1e-14,gradient) + 2019-04-05T14:18:18: Index ....................: -3.21741632072522 +/- 0.0232805723886611 [-5,5] (free,scale=1,gradient) + 2019-04-05T14:18:18: PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient) + 2019-04-05T14:18:18: Number of temporal par's ..: 1 + 2019-04-05T14:18:18: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2019-04-05T14:18:18: === GCTAModelAeffBackground === + 2019-04-05T14:18:18: Name ......................: bkg_23559 + 2019-04-05T14:18:18: Instruments ...............: HESS + 2019-04-05T14:18:18: Instrument scale factors ..: unity + 2019-04-05T14:18:18: Observation identifiers ...: 23559 + 2019-04-05T14:18:18: Model type ................: "PowerLaw" * "Constant" + 2019-04-05T14:18:18: Number of parameters ......: 4 + 2019-04-05T14:18:18: Number of spectral par's ..: 3 + 2019-04-05T14:18:18: Prefactor ................: 2.26357632227583e-13 +/- 3.56983854227473e-15 [1e-16,1e-12] ph/cm2/s/MeV (free,scale=1e-14,gradient) + 2019-04-05T14:18:18: Index ....................: -3.25328282552122 +/- 0.0239266148104449 [-5,5] (free,scale=1,gradient) + 2019-04-05T14:18:18: PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient) + 2019-04-05T14:18:18: Number of temporal par's ..: 1 + 2019-04-05T14:18:18: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2019-04-05T14:18:18: === GCTAModelAeffBackground === + 2019-04-05T14:18:18: Name ......................: bkg_23592 + 2019-04-05T14:18:18: Instruments ...............: HESS + 2019-04-05T14:18:18: Instrument scale factors ..: unity + 2019-04-05T14:18:18: Observation identifiers ...: 23592 + 2019-04-05T14:18:18: Model type ................: "PowerLaw" * "Constant" + 2019-04-05T14:18:18: Number of parameters ......: 4 + 2019-04-05T14:18:18: Number of spectral par's ..: 3 + 2019-04-05T14:18:18: Prefactor ................: 2.87721724170504e-13 +/- 4.31204359584396e-15 [1e-16,1e-12] ph/cm2/s/MeV (free,scale=1e-14,gradient) + 2019-04-05T14:18:18: Index ....................: -3.25393054001656 +/- 0.0227802144778143 [-5,5] (free,scale=1,gradient) + 2019-04-05T14:18:18: PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient) + 2019-04-05T14:18:18: Number of temporal par's ..: 1 + 2019-04-05T14:18:18: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2019-04-05T14:18:18: === GModelSky === + 2019-04-05T14:18:18: Name ......................: Crab + 2019-04-05T14:18:18: Instruments ...............: all + 2019-04-05T14:18:18: Instrument scale factors ..: unity + 2019-04-05T14:18:18: Observation identifiers ...: all + 2019-04-05T14:18:18: Model type ................: PointSource + 2019-04-05T14:18:18: Model components ..........: "PointSource" * "PowerLaw" * "Constant" + 2019-04-05T14:18:18: Number of parameters ......: 6 + 2019-04-05T14:18:18: Number of spatial par's ...: 2 + 2019-04-05T14:18:18: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) + 2019-04-05T14:18:18: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) + 2019-04-05T14:18:18: Number of spectral par's ..: 3 + 2019-04-05T14:18:18: Prefactor ................: 1.205595290581e-15 +/- 1.02399444013702e-16 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) + 2019-04-05T14:18:18: Index ....................: -2.64310065863216 +/- 0.0532926238637032 [-0,-5] (free,scale=-1,gradient) + 2019-04-05T14:18:18: PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient) + 2019-04-05T14:18:18: Number of temporal par's ..: 1 + 2019-04-05T14:18:18: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2019-04-05T14:18:18: === GCTAModelIrfBackground === + 2019-04-05T14:18:18: Name ......................: CTABackgroundModel + 2019-04-05T14:18:18: Instruments ...............: CTA + 2019-04-05T14:18:18: Instrument scale factors ..: unity + 2019-04-05T14:18:18: Observation identifiers ...: all + 2019-04-05T14:18:18: Model type ................: "PowerLaw" * "Constant" + 2019-04-05T14:18:18: Number of parameters ......: 4 + 2019-04-05T14:18:18: Number of spectral par's ..: 3 + 2019-04-05T14:18:18: Prefactor ................: 1 +/- 0 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) + 2019-04-05T14:18:18: Index ....................: 0 +/- 0 [-5,5] (free,scale=1,gradient) + 2019-04-05T14:18:18: PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient) + 2019-04-05T14:18:18: Number of temporal par's ..: 1 + 2019-04-05T14:18:18: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) .. note:: + Using ``inmodel=$CTOOLS/share/models/crab.xml`` will also append the CTA + background model to the ``models.xml`` file, yet since this model only + applies to CTA data it is ignored in the fit. This is signalled by - By default, :ref:`csfindobs` only selects data of highest quality (i.e. QUALITY=0). - You may overwrite this default by specifying the hidden parameter - ``min_qual``. For example, ``min_qual=1`` selects all data with a - looser quality criteria. + .. code-block:: none + 2019-04-05T14:18:17: Parameter "Prefactor" has zero curvature. Fix parameter. + 2019-04-05T14:18:17: Parameter "Index" has zero curvature. Fix parameter. -Create an observation list --------------------------- - -As next step, you must convert the runlist ASCII file into an observation -definition XML file. The observation definition XML file will contain the -file names of all files that are needed for the analysis. You do this -conversion with the :ref:`csiactobs` script: - -.. code-block:: bash - - $ csiactobs - Data storage name [] fits-data-name - Input runlist file [runlist.lis] - Number of free parameters per background model [1] - Output model definition XML file [bgmodels.xml] - Output observation definition XML file [obs.xml] - -The :ref:`csiactobs` script will create two output files: the observations -definition XML file ``obs.xml`` and an output model definition XML file -``bgmodels.xml``. To generate ``obs.xml``, :ref:`csiactobs` has used the -IACT data storage and extracted the relevant file names. ``bgmodels.xml`` is -a file that is used for background modeling, where each observation will have -its own independent background model. In the example above, you have set the -number of free parameters per background model to one, hence the normalisation -of the background model for each observation will be a free parameter that is -later adjusted by the maximum likelihood fit. - -There are some further hidden parameters to steer the start parameters for the -effective area and the Gaussian background models. You may have a look at -:ref:`csiactobs` to see the full list of available parameters. - -You may also note that the :ref:`csiactobs` script has created a ``csiactobs.log`` -file that logs the complete energy range of the observations that have been -selected. These values may be important for later, in particular if you -intend to do a stacked analysis. - -In case that you have already a model definition XML file that describes a model -of the celestial source distribution (a.k.a. a sky model), you may provide this -sky model to :ref:`csiactobs` using the hidden ``inmodel`` parameter. The -output model definition XML file will then contain the sky model and the -background model, and you can use the model directly in a maximum likelihood -fit. - -.. code-block:: bash - - $ csiactobs inmodel="mymodel.xml" - -Alternatively, you can merged models using the :ref:`csmodelmerge` script: - -.. code-block:: bash - - $ csmodelmerge - Input model definition XML files [mymodel.xml bgmodels.xml] - Output model definition XML file [combined_model.xml] - -.. note:: - - The number of files that can be merged with :ref:`csmodelmerge` is not - limited to two. The list of input file names may be either separated by - whitespace or semi-colons, can be specified using wildcards, or can be given - in an ASCII file (see :ref:`csmodelmerge`). + during the fit, and the corresponding model parameter have an error of + zero, although they are signalled as free. .. warning:: - - Each model component in the input model definition XML files that are provided - to :ref:`csmodelmerge` must have a different and unique name. Merging for - example the same XML model twice will raise an exception. - -A list of available sky models can be found -`here `_. -If you are not familiar with creating sky models you should read the section -about :ref:`modelling CTA data `. - - -Example XML files ------------------ -To get familiar with the XML syntax and format, example files for an observation -definition XML file and a model definition XML file are shown in the following. - -Observation definition XML file -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. code-block:: xml - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -Model definition XML file -^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. code-block:: xml - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -.. note:: - - It is important to ensure background models are properly linked to their - respective observation. Therefore it is required to keep the attributes - ``instrument`` and ``id`` the same for the observation and the corresponding - background model. The tool :ref:`csiactobs` assures this automatically. - - -Run ctselect ------------- -To prepare the data for analysis, cuts have to be applied to the event data. -The selection is performed by :ref:`ctselect`. -This tool writes out selected event lists into the local directory. If the -observation definition XML file contains several runs, it is recommended -to first create a separate folder and specify this folder in the hidden -``prefix`` argument. - -.. code-block:: bash - - $ mkdir selected - -.. code-block:: bash - - $ ctselect usethres=DEFAULT usepnt=yes prefix=selected/selected_ - Input event list or observation definition file [events.fits] obs.xml - Lower energy limit (TeV) [0.1] - Upper energy limit (TeV) [100.0] - Radius of ROI (degrees) (0-180) [3.0] 2.5 - Start time (CTA MET in seconds) [0.0] - End time (CTA MET in seconds) [0.0] - Output event list or observation definition file [selected_events.fits] selected_obs.xml - -For IACT analysis, it is recommended to use the hidden parameter -``usethres="DEFAULT"``. This instructs :ref:`ctselect` to extract the safe -energy range from the instrument response functions and apply them to the -data. This safe energy range is thus superior to the energy limit passed via -the user parameters. In addition, to analyse the complete field of view, the -parameter ``usepnt=yes`` uses, for each observation, the pointing position as -centre for the selection radius. The radius parameter is dependent on the -intrument, for an instrument with a 5° field of view, a radius of 2.5° -seems reasonable. The time selection is not applied in the above example; -specifying 0 as start and end time skips the time selection. For time-resolved -analysis, it is important to know the MET time that is required to extract. -The result of the selection step is written into the observation XML file -``selected_obs.xml``, which now contains references to the new selected event -FITS files. - - -Unbinned analysis ------------------ -Once the data is selected, the easiest way to analyse is an unbinned analysis. -Note that the input model definition XML file must now contain the background -models and source components to describe the field of view. - -.. code-block:: bash - - $ ctlike - Event list, counts cube or observation definition file [events.fits] selected_obs.xml - Source model [$CTOOLS/share/models/crab.xml] crab_models.xml - Source model output file [crab_results.xml] - -The result of the fit was stored in ``crab_results.xml``. Note that fitted -parameters, ``Prefactors`` in particular, typically use MeV as energy unit. -To monitor the progress of the fit on the screen, one can simply run with the -option ``debug=yes``. Alternatively, the logfile ``ctlike.log`` can be inspected -after the fit. - -On default, energy dispersion is not considered in the fit. To switch on the -usage of the energy migration matrix, the hidden parameter ``edisp=yes`` can -be provided. Note that this will cause a significant reduction of the -computing speed. - - -Stacked analysis ----------------- -In a stacked analysis the events of all observations are stacked into a single -counts cube. This analysis mode is much faster than unbinned analysis when -having a large dataset (e.g. > 100 hours). For this type of analysis, some -intermediate data products have to be produced. The products are a binned counts -cube, an exposure cube, a point spead function cube, optionally an energy -dispersion cube, and a background cube. - -Bin events -^^^^^^^^^^ - -.. code-block:: bash - - $ ctbin - Event list or observation definition file [events.fits] selected_obs.xml - 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] - Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] - Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] - Image scale (in degrees/pixel) [0.02] - Size of the X axis in pixels [200] - Size of the Y axis in pixels [200] - Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] - Start value for first energy bin in TeV [0.1] 0.5 - Stop value for last energy bin in TeV [100.0] 50 - Number of energy bins [20] - Output counts cube [cntcube.fits] - - -Create exposure cube -^^^^^^^^^^^^^^^^^^^^ -After binning the events into a three-dimensional cube, an exposure cube has to -be computed. The exposure is defined as the effective area times the dead-time -corrected observation time. Each observation from the input container gets -stacked in the resulting cube. The exposure is stored in units of :math:`cm^2 s`. -The exposure cube does not have to contain the same binning as the event cube -but for simplicity, the event cube can be passed to adopt the binning -parameters. Note, however, that the exposure cube is defined in true sky -coordinates and energy while the counts cube is defined in reconstructed sky -coordinates and energy. Consequently, the sky area and energy range covered by -the exposure cube should be slightly larger than that of the counts cube to -accommodate for spill over of events due to the point spread function and energy -dispersion. By computing the exposure cube on the same grid as the counts cube, -the spill over of events from sources at the edge of cube will not be handled -correctly. In the example below, however, no source at the edge of the field of -view is present. Therefore, for simplicity, the count cube is used as input cube -to extract the binning. - -This task of computing the exposure cube is is performed by :ref:`ctexpcube`. - -.. code-block:: bash - - $ ctexpcube - Event list or observation definition file [NONE] selected_obs.xml - Input counts cube file to extract exposure cube definition [NONE] cntcube.fits - Output exposure cube file [expcube.fits] - -Alternatively, the exposure cube can be created with different binning than the -counts cube: - -.. code-block:: bash - - $ ctexpcube - Input event list or observation definition XML file [NONE] selected_obs.xml - Input counts cube file to extract exposure cube definition [NONE] - 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] - Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] - Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] - Image scale (in degrees/pixel) [0.02] 0.04 - Size of the X axis in pixels [200] 100 - Size of the Y axis in pixels [200] 100 - Lower energy limit (TeV) [0.5] - Upper energy limit (TeV) [50.0] - Number of energy bins [20] 30 - Output exposure cube file [expcube.fits] - -Create point spread function cube -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -As a next step for the stacked analysis, a cube containing the point spread -function (PSF) must be computed. Since the PSF cannot be stored by one single -parameter, the PSF cube computed by :ref:`ctpsfcube` has a fourth dimension. -In each bin of the cube, the PSF is stored as a function of offset from source. -The granularity of the PSF histogram is determined by the hidden parameter -``anumbins`` (default: 200). Therefore, when passing the event cube to adopt -the sky binning for the PSF cube, the resulting FITS file can become quite large -due to the fourth dimension. Usually in IACT analysis, the PSF doesn't change -too dramatically across the field of view. Therefore the user can think about -reducing the spatial binning of the PSF cube: - -.. code-block:: bash - - $ ctpsfcube - Input event list or observation definition XML file [NONE] selected_obs.xml - Input counts cube file to extract PSF cube definition [NONE] - 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] - Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] - Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] - Image scale (in degrees/pixel) [1.0] 0.2 - Size of the X axis in pixels [10] 20 - Size of the Y axis in pixels [10] 20 - Lower energy limit (TeV) [0.1] 0.5 - Upper energy limit (TeV) [100.0] 50 - Number of energy bins [20] - Output PSF cube file [psfcube.fits] - -Depending on the required PSF precision, one could reduce the number of offset -bins via the hidden parameter ``anumbins``: - -.. code-block:: bash - - $ ctpsfcube anumbins=100 - -Create energy dispersion cube (optional) -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -In case energy dispersion should be considered used in the analysis, the tool -:ref:`ctedispcube` computes the energy migration response for the given -observations. Analoguous to the PSF cube, the energy dispersion cube is stored -in a four-dimensional map. In each bin of the cube, the migration value -``E_reco/E_true`` is stored. The granularity of the migration histogram can be -steered via the hidden parameter ``migrabins`` (default: 100). In addition, the -maximum migration value can be set via the hidden parameter ``migramax`` -(default: 2.0). Similar to the PSF cube, the energy dispersion cube FITS file -can become quite large. Therefore, since the energy dispersion shouldn't vary -too much across the sky region, one could reduce the spatial binning: - -.. code-block:: bash - - $ ctedispcube - Input event list or observation definition XML file [NONE] selected_obs.xml - Input counts cube file to extract energy dispersion cube definition [NONE] - 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] - Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] - Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] - Image scale (in degrees/pixel) [1.0] 0.2 - Size of the X axis in pixels [10] 20 - Size of the Y axis in pixels [10] 20 - Lower energy limit (TeV) [0.1] 0.5 - Upper energy limit (TeV) [100.0] 50 - Number of energy bins [20] - Output energy dispersion cube file [edispcube.fits] - -Create background cube -^^^^^^^^^^^^^^^^^^^^^^ -Last but not least, for a stacked IACT analysis a cube containing the background -rate in sky coordinates and reconstructed energy has to be computed. This task -is performed by :ref:`ctbkgcube`. The binning here can also differ from the counts -cube. For simplicity, however, the example below uses the counts cube to adopt -the binning. - -.. code-block:: bash - - $ ctbkgcube debug=yes - Input event list or observation definition XML file [NONE] selected_obs.xml - Input counts cube file to extract background cube definition [NONE] cntcube.fits - Input model XML file [NONE] crab_models.xml - Output background cube file [bkgcube.fits] - Output model XML file [NONE] binned_models.xml - -Note that this tool also requires the parameters of an input and output model. -In the model XML file that came out of :ref:`csiactobs`, one background model -per observation is included. This models get merged and averaged in the background -sky cube. In the output model the background models per obseervation will be -removed. Instead, a global background model for the newly created background -cube is included. Sky models present in the input model XML file will also be -included in the new XML file, which subsequently can be used for stacked -:ref:`ctlike`. - -Example for stacked model XML file -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The output model of :ref:`ctbkgcube` looks the following: - -.. code-block:: xml - - - - - - - - - - - - - - - - - - - - - - - -The background model with ``type=CTACubeBackground`` is used to scale the -background cube stored in the FITS file created by :ref:`ctbkgcube`. - -Run ctlike -^^^^^^^^^^ -Having all the intermediate data products ready, a stacked analysis can be -conducted using :ref:`ctlike`. - -.. code-block:: bash - - $ ctlike - Input event list, counts cube or observation definition XML file [events.fits] cntcube.fits - Input exposure cube file (only needed for stacked analysis) [NONE] expcube.fits - Input PSF cube file (only needed for stacked analysis) [NONE] psfcube.fits - Input background cube file (only needed for stacked analysis) [NONE] bkgcube.fits - Input model XML file [binned_models.xml] - Output model XML file [binned_results.xml] - -Note that when passing an counts cube to :ref:`ctlike`, the tool behaves -differently than in unbinned mode. It queries directly for the additional -ingredients for the stacked analysis. It is important to pass the background -model generated by :ref:`ctbkgcube` here to ensure the proper modelling of -the background in the fit. - -To consider also the energy dispersion, you have to execute the :ref:`ctlike` -tool with the hidden parameter ``edisp=yes``. The tool will then also query for -the energy dispersion cube: - -.. code-block:: bash - - $ ctlike debug=yes - Input event list, counts cube or observation definition XML file [events.fits] cntcube.fits - Input exposure cube file (only needed for stacked analysis) [NONE] expcube.fits - Input PSF cube file (only needed for stacked analysis) [NONE] psfcube.fits - Input energy dispersion cube file (only needed for stacked analysis) [NONE] edispcube.fits - Input background cube file (only needed for stacked analysis) [NONE] bkgcube.fits - Input model XML file [binned_models.xml] - Output model XML file [binned_results.xml] + The effective area background model generated by :ref:`csiactobs` is very + simplistic, and should not be used for a serious data analysis of H.E.S.S. + data. Please follow the procedure :ref:`hess_dr1` for the generation of + a more reliable background model. diff --git a/doc/source/users/tutorials/howto/iact/copy.rst b/doc/source/users/tutorials/howto/iact/copy.rst index 5ac7365..3b3c57c 100755 --- a/doc/source/users/tutorials/howto/iact/copy.rst +++ b/doc/source/users/tutorials/howto/iact/copy.rst @@ -3,66 +3,162 @@ Copying IACT data ================== + .. admonition:: What you will learn + + You will learn how copy data from an IACT data store to your disk. + Before you start ---------------- -Currently the way how IACT data is stored is quite prototypical. The ctools can work -with IACT data that is stored following the `specifications here `_. -The script :ref:`csiactcopy` can be used to transfer these data from one system to another. -IACT data is the property of collaborations (such as HESS, MAGIC or VERITAS) and not publicly available. If you are member of one of these -collaborations, you should have access to a central place where the FITS data is stored and maintained. - -Mount remote file system ------------------------- -Note that :ref:`csiactcopy` does not include functionality to establish a conection to a remote file system. -For the moment, it is capable of copying data from one location to another. Therefore it is necessary to -mount the remote file system where the data is located on the local machine. This can be achieved using e.g. -`sshfs `_: + +Currently the way how IACT data is stored is quite prototypical. +The ctools can work with IACT data that is stored following the specifications +provided in the +`Data formats for gamma-ray astronomy `_. +The script :ref:`csiactcopy` can be used to transfer IACT data from one system +to another. +IACT data is the property of collaborations (such as HESS, MAGIC or VERITAS) +and not publicly available, except for the recent +`public release of some H.E.S.S. data `_. +In this tutorial we use these public data. +After downloading the data and uncompression the tarball you should have a +folder named ``hess_dl3_dr1`` in your current workingcdirectory with the +following structure: .. code-block:: bash - $ sshfs user@remote.server.com:/path/to/remote/fits/data/ /path/to/mountpoint/ + hess_dl3_dr1/ + hess_dl3_dr1/data + hess_dl3_dr1/hdu-index.fits.gz + hess_dl3_dr1/obs-index.fits.gz + hess_dl3_dr1/README.txt + +You have to add the file ``master.json`` to this folder with the following +content: + +.. code-block:: json + + { + "datasets": [ + { + "name": "hess_dl3_dr1", + "hduindx": "hdu-index.fits.gz", + "obsindx": "obs-index.fits.gz", + "Total observations": 106, + "Total time": "48.6 hours", + "Total events": 662826, + "Disk usage": "43 MB", + "Date created": "2018-09-10" + } + ] + } + +.. note:: + If you are member of an IACT collaboration and if the collaboration has + setup a datastore in the + `compliant format `_ + you can also establish a remote connection to this datastore using for + example + + .. code-block:: bash + + $ sshfs user@remote.server.com:/path/to/remote/fits/data/ /project-data/hess + + The success of this operation can be verified by checking if there is a file + called ``master.json`` in the mounted directory. + + +Check what data are available +----------------------------- + +Before you start an analysis, it is important that you know what kind of data +are available for analysis. +In the `current scheme `_ +every dataset is represented by a unique name, in our case ``hess_dl3_dr1``. +You have to pass this string to some set of scripts to find and use the data +for analysis. + +To check what datasets are available you should run the :ref:`csiactdata` +script: + +.. code-block:: bash + + $ csiactdata debug=yes + Path were data is located [] /project-data/hess/hess_dl3_dr1 + +This script will list the names of available datasets. Here at excerpt of the +log file ``csiactdata.log``: + +.. code-block:: none + + 2019-04-05T13:10:35: +======================+ + 2019-04-05T13:10:35: | Data storage entries | + 2019-04-05T13:10:35: +======================+ + 2019-04-05T13:10:35: Master index file .........: /project-data/hess/hess_dl3_dr1/master.json + 2019-04-05T13:10:35: + 2019-04-05T13:10:35: +------------------------+ + 2019-04-05T13:10:35: | Available data configs | + 2019-04-05T13:10:35: +------------------------+ + 2019-04-05T13:10:35: === hess_dl3_dr1 === + 2019-04-05T13:10:35: Name ......................: hess_dl3_dr1 + 2019-04-05T13:10:35: Observation index .........: obs-index.fits.gz + 2019-04-05T13:10:35: HDU index .................: hdu-index.fits.gz + 2019-04-05T13:10:35: Date created ..............: 2018-09-10 + 2019-04-05T13:10:35: Total events ..............: 662826 + 2019-04-05T13:10:35: Total time ................: 48.6 hours + 2019-04-05T13:10:35: Disk usage ................: 43 MB + 2019-04-05T13:10:35: Total observations ........: 105 -The success of this operation can be verified by checking if there is a file called ``master.json`` -in the mounted directory. Copy data --------- -After mounting the file system to a folder of your choice (``/path/to/mountpoint/`` in this case), the copying of the data can be initiated. -To monitor the progress on the screen, the hidden parameter ``debug=yes`` can be used. To increase the chattiness of the output, -the hidden parameter ``chatter`` can e.g. be set to 3. -.. code-block:: bash - - $ csiactcopy debug=yes chatter=3 - Location of remote master file [/path/to/mountpoint/master.json] - Name of FITS production to download [iact-fits] - Destination path of FITS data [/path/to/local/fits/data/] - -The name of the FITS production must be available in the remote -``master.json`` file. In case the names of the available productions are unkown -there are two options to determine an appropriate name: - -* Run :ref:`csiactdata` and set as ``datapath`` the mountpoint -* Try any name, :ref:`csiactcopy` will throw a ``RuntimeError`` with the available production names if the input name is not found - -Copy only a subset of the data ------------------------------- -If the user is interested only in a small subset of the remote FITS production, -it is also possible to download only parts of the data. The user simply has to provide an -ASCII file containing a list of observation IDs to be downloaded (one ID per line). The hidden -parameter ``runlist`` then needs to be specified. +You can now copy the data from the data store to your disk. In this example we +assume that the H.E.S.S. datastore is at ``/project-data/hess/hess_dl3_dr1``. +To do the copy you should type .. code-block:: bash - - $ csiactcopy runlist=myrunlist.txt debug=yes - + + $ csiactcopy + Location of remote master file [/path/to/mountpoint/master.json] /project-data/hess/hess_dl3_dr1/master.json + Name of FITS production to download [iact-fits] hess_dl3_dr1 + Destination path of FITS data [/path/to/local/fits/data/] data + +.. note:: + To monitor the progress on the screen you should set the hidden parameter + ``debug=yes``. You may also increase the chattiness of the output by + increasing the hidden parameter ``chatter`` from its default value of 2 + to for example 3. + +.. note:: + If you are only interested in a subset of the dataset you can also download + only that subset. For that purpose you have to specify an ASCII file + containing a list of observation identifiers to be downloaded, e.g. + + .. code-block:: none + + 20136 + 20137 + 20151 + + You can then run :ref:`csiactcopy` by specifying the ASCII filename + for the hidden ``runlist`` parameter: + + .. code-block:: bash + + $ csiactcopy runlist=runlist.txt + Location of remote master file [/project-data/hess/hess_dl3_dr1/master.json] + Name of FITS production to download [hess_dl3_dr1] + Destination path of FITS data [data] + + Troubleshooting --------------- -In case a download was interupted (e.g. the remote file system was disconnected), the script can simply be rerun. -The index files which list the files available on the system are refreshed after the copying process. -The hidden boolean parameter ``clobber`` (default=yes) specifies if the user wishes that local content should be overwritten with -remote content. In turn, specifying ``clobber=no`` can speed up the process since remote content does not need to be copied -if it is already available on the user machine. - -To verify the copying process was successful, it is recommended to run :ref:`csiactdata` after the process in order to list available productions -on the local machine. The output should now include the copied FITS production name. +In case a download was interupted (e.g. the remote file system was +disconnected), the script can simply be rerun. +The index files, which list the files available on the system, are refreshed +after the copying process. +The hidden boolean parameter ``clobber`` (by default ``yes``) specifies if the +user wishes that local content should be overwritten with remote content. +In turn, specifying ``clobber=no`` can speed up the process since remote +content does not need to be copied if it is already available on the user +machine. diff --git a/doc/source/users/tutorials/howto/iact/high_level.rst b/doc/source/users/tutorials/howto/iact/high_level.rst deleted file mode 100755 index 4550691..0000000 --- a/doc/source/users/tutorials/howto/iact/high_level.rst +++ /dev/null @@ -1,315 +0,0 @@ -.. _sec_iact_high_level: - -High level analysis tools for IACT data -======================================= - -This section describes and explains how to create high level output. The -following scripts and tools are explained: - -* :ref:`csobsinfo` (inspect and observation XML files) -* :ref:`csmodelinfo` (inspect model XML files) -* :ref:`ctulimit` (compute upper limits) -* :ref:`cterror` (compute asymmetric errors) -* :ref:`csspec` (compute spectral points) -* :ref:`cslightcrv` (compute a light curve) -* :ref:`csresmap` (compute a residual map) -* :ref:`cttsmap` (compute a test statistics map) - -Inspecting observation definition files ---------------------------------------- -After having created an observation definition file, the user might want to -know a summary of its content. Important quantities usually are the number of -observations, the total livetime, or the mean zenith angle. For this purpose, -:ref:`csobsinfo` inspects an observation XML container for its content. - -.. code-block:: bash - - $ csobsinfo - Event list, counts cube, or observation definition file [obs.xml] selected_obs.xml - -This dumps general information of the observation container in a logfile ``csobsinfo.log``. -For a direct output on the screen, pass the hidden parameter ``debug=yes`` on the command line: - -.. code-block:: bash - - $ csobsinfo debug=yes - Event list, counts cube, or observation definition file [obs.xml] selected_obs.xml - ... - 2016-02-19T15:43:08: +=========+ - 2016-02-19T15:43:08: | Summary | - 2016-02-19T15:43:08: +=========+ - 2016-02-19T15:43:08: Unbinned observations .....: 4 - 2016-02-19T15:43:08: Number of events ..........: 6085 - 2016-02-19T15:43:08: Binned observations .......: 0 - 2016-02-19T15:43:08: Number of bins ............: 0 - 2016-02-19T15:43:08: - 2016-02-19T15:43:08: === Pointings === - 2016-02-19T15:43:08: Mean zenith angle .........: 46.89 - 2016-02-19T15:43:08: Mean azimuth angle ........: 13.63 - 2016-02-19T15:43:08: - 2016-02-19T15:43:08: === Energy range === - 2016-02-19T15:43:08: Emin ......................: 461.48 GeV - 2016-02-19T15:43:08: Emax ......................: 37.8832 TeV - 2016-02-19T15:43:08: - 2016-02-19T15:43:08: === Time range === - 2016-02-19T15:43:08: Start [MJD] ...............: 53343.9212153 - 2016-02-19T15:43:08: Stop [MJD] ................: 53347.9321065 - 2016-02-19T15:43:08: Total ontime ..............: 6877.00 s = 114.62 min = 1.91 h - 2016-02-19T15:43:08: Total livetime ............: 6427.52 s = 107.13 min = 1.79 h - ... - -Further options are described on the reference page :ref:`csobsinfo`. - -Inspecting model XML files --------------------------- -Similar to the observation containers, model XML file can be inspected, too. -The tool :ref:`csmodelinfo` gives a summary of the model container. In particular, -the number of free parameters, number of sky and background models and also number of -parameters that are at their limits might be of interest. - -.. code-block:: bash - - $ csmodelinfo debug=yes - Input model definiton file XML file [$CTOOLS/share/models/crab.xml] - Output DS9 region file [ds9.reg] - 2016-02-19T15:50:40: +=========+ - 2016-02-19T15:50:40: | Summary | - 2016-02-19T15:50:40: +=========+ - 2016-02-19T15:50:40: === Instrument specific models === - 2016-02-19T15:50:40: All .......................: 1 - 2016-02-19T15:50:40: CTA .......................: 1 - 2016-02-19T15:50:40: === Model types === - 2016-02-19T15:50:40: PointSource ...............: 1 - 2016-02-19T15:50:40: CTAIrfBackground ..........: 1 - 2016-02-19T15:50:40: - 2016-02-19T15:50:40: +=======================+ - 2016-02-19T15:50:40: | Parameter information | - 2016-02-19T15:50:40: +=======================+ - 2016-02-19T15:50:40: All parameters ............: 10 - 2016-02-19T15:50:40: Fixed parameters ..........: 6 - 2016-02-19T15:50:40: Free parameters (total) ...: 4 - 2016-02-19T15:50:40: Free background parameters : 2 - 2016-02-19T15:50:40: Free source parameters ....: 2 - 2016-02-19T15:50:40: Free spectral parameters ..: 2 - 2016-02-19T15:50:40: Parameters at limit .......: 0 - -Similar to :ref:`csobsinfo`, the script also generates a ds9 region file -including the positions of all the sky model components. - -Compute upper limit -------------------- -Very often in gamma-ray astronomy sources are at the verge of detection or even not detectable. -In such cases, it is useful to derive an upper limit using :ref:`ctulimit`. - -.. code-block:: bash - - $ ctulimit - Input event list, counts cube or observation definition XML file [selected_obs.xml] - Input model XML file [crab_models.xml] - Source of interest [Crab] - -The upper limit will be stored in the log file. To get the limit printed on screen, use the hidden parameter -``debug=yes``. - -Compute asymmetric errors -------------------------- -When an analysis approaches its final state, it makes sense to have asymmetric errors on the parameters of the -source of interest. For this purpose, the tool :ref:`cterror` can be used: - -.. code-block:: bash - - $ cterror debug=yes - Input event list, counts cube or observation definition XML file [selected_obs.xml] - Input model XML file [crab_models.xml] - Source of interest [Crab] - Output model XML file [cterror_results.xml] - -The output model does not contain asymmetric errors yet. The positive and negative uncertainties can be read from the -logfile (or from screen if ``debug=yes`` was specified). - -Compute spectral points ------------------------ -A very common task in astronomy is to compute spectral data points. To determine a spectral data point, a small energy -range is considered and the model prefactor and its uncertainty is evaluated. The tool :ref:`csspec` works on both, binned -and unbinned data. The hidden parameter ``edisp=yes`` can be specified in both cases to consider the energy migration -matrix in the fit. - -.. note:: - - The resulting spectral points are provided as a function of reconstructed energy. - -Unbinned -^^^^^^^^ - -.. code-block:: bash - - $ csspec debug=yes - Input event list, counts cube, or observation definition XML file [events.fits] selected_obs.xml - Input model XML file [$CTOOLS/share/models/crab.xml] crab_models.xml - Source name [Crab] - Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] - Lower energy limit for spectral points (TeV) [0.1] 0.5 - Upper energy limit for spectral points (TeV) [100.0] 50.0 - Number of spectral points (1-10000) [20] 10 - Output spectrum file [spectrum.fits] - -Binned -^^^^^^ - -.. code-block:: bash - - $ csspec debug=yes - Input event list, counts cube, or observation definition XML file [cntcube.fits] - Input exposure cube file (only needed for stacked analysis) [expcube.fits] - Input PSF cube file (only needed for stacked analysis) [psfcube.fits] - Input background cube file (only needed for stacked analysis) [bkgcube.fits] - Input model XML file [binned_models.xml] - Source name [Crab] - Number of spectral points (1-10000) [10] - Lower energy limit for spectral points (TeV) [0.5] - Upper energy limit for spectral points (TeV) [50.0] - Output spectrum file [spectrum.fits] - -Plot spectral points -^^^^^^^^^^^^^^^^^^^^ -Instead of supporting a plotting library, simple example scripts are available to visualise data products. -Have a look at ``$CTOOLS/examples/show_spectrum.py`` how a spectrum can be plotted using ``matplotlib``. -To have first glance at the above computed spectrum one can use this script in the following way: - -.. code-block:: bash - - $ python $CTOOLS/examples/show_spectrum.py spectrum.fits - - -Compute light curves --------------------- -Euqally common in astronomy are light curves from a time-series analysis. The tool :ref:`cslightcrv` takes as input only -an unbinned observation container (or a single event list). Counts cube are not possible as inout since the time informstion -is lost during the binning procedure. Nevertheless, in case of large time bins, it is possible to require a binned analysis -in each time bin. The tool will accordingly slice the event list and create the data products for the stacked analysis in each -time span. - -.. code-block:: bash - - $ cslightcrv - Input event list, counts cube, or observation definition XML file [selected_obs.xml] - Input model XML file [crab_models.xml] - Source name [Crab] - Algorithm for defining time bins (FILE|LIN|GTI) [GTI] - Number of energy bins per light curve bin (0=unbinned) [0] - Lower energy limit of events (TeV) [0.1] 0.5 - Upper energy limit of events (TeV) [100.0] 50 - Output light curve file [lightcurve.fits] - -Compute a residual map ----------------------- -A frequent means to visually inspect the fitted model with respect to the input data is to create residual maps. -Using :ref:`csresmap` such a map can easily be computed. The tool internally bins the data according to user parameters. -(if data is not already provided in a binned state). Taking into account the corresponding IRFs, -:reF:`csresmap` computes a model map (running :ref:`ctmodel`). Subsequently, there are several choices how data -and model should be compared. There are three options: - -* ``SUB``: the subtraction of the model from the counts. The resulting map will display differences in absolute counts. -* ``SUBDIV``: the subtraction and division by the model. The resulting map will display relative differences of the data with respect to the model. -* ``SUBDIVSQRT``: the subtraction and division by the square root of the model. The resulting map will display an approximation of a residual significance map. In case of sufficient count statistic per bin, the bin value represents the significance. - -Example unbinned observation container -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. code-block:: bash - - $ csresmap - Input event list, counts cube, or observation definition XML file [selected_obs.xml] - Input model XML file [crab_results.xml] - 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] - Coordinate System (CEL|GAL) [CEL] - Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] - Size of the X axis in pixels [100] - Size of the Y axis in pixels [100] - Pixel size (deg/pixel) [0.02] - Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUBDIV] - Output residual map file [resmap.fits] - -Example binned/stacked observation container -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The tool behaves differently if the input observation parameter is a count cube. It will query for the additional binned reponse components. - -.. code-block:: bash - - $ csresmap - Input event list, counts cube, or observation definition XML file [cntcube.fits] - Input model cube file (generated with ctmodel) [NONE] - Input exposure cube file (only needed for stacked analysis) [expcube.fits] - Input PSF cube file (only needed for stacked analysis) [psfcube.fits] - Input background cube file (only needed for stacked analysis) [bkgcube.fits] - Input model XML file [binned_results.xml] - Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUBDIV] - Output residual map file [resmap.fits] - -Note that the tool also queries for an input model cube file (which we set to ``NONE`` here). This is very convenient in case the model cube has already been precomputed -using :ref:`ctmodel`. - -Example with input model cube -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -In case a model cube can be provided, :ref:`csresmap` simply collapses both cube into one skymap and applies the specified algorithm. - -.. code-block:: bash - - $ csresmap - Input event list, counts cube, or observation definition XML file [cntcube.fits] - Input model cube file (generated with ctmodel) [modcube.fits] - Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUB] - Output residual map file [resmap.fits] - - -Compute a test statistics (TS) map ----------------------------------- -The test statistics quantity is a more precise measurement of statistic than the ``SUBDIVSQRT`` option in :ref:`csresmap`. -To get familiar with the concept of this procedure, :ref:`read more about TS maps `. - -.. code-block:: bash - - $ cttsmap - Input event list, counts cube or observation definition XML file [selected_obs.xml] - Input model XML file [crab_results.xml] - Test source name [Crab] - 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] - Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] - Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] - Image scale (in degrees/pixel) [0.05] - Size of the X axis in pixels [20] - Size of the Y axis in pixels [20] - Output Test Statistic map file [tsmap.fits] - -The output file will contain a sky map holding the TS map. In complementary extensions, for each free parameter of the source model, a map displaying the parameter -value is added to the output file. For instance a map of the prefactor can be visualised with `DS9 `_ as follows: - -.. code-block:: bash - - $ ds9 tsmap.fits[Prefactor] - -Since :ref:`cttsmap` has to run a parameter optimisation in each skymap bin, it is very time comsuming to -compute a fine-granulated TS sky map. Read here how to split up the :ref:`computation into several jobs `. -Similar as many other tools :ref:`cttsmap` can also work on binned observation input: - -.. code-block:: bash - - $ cttsmap debug=yes - Input event list, counts cube or observation definition XML file [cntcube.fits] - Input exposure cube file (only needed for stacked analysis) [expcube.fits] - Input PSF cube file (only needed for stacked analysis) [psfcube.fits] - Input background cube file (only needed for stacked analysis) [bkgcube.fits] - Input model XML file [binned_results.xml] - Test source name [Crab] - 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] - Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] - Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] - Image scale (in degrees/pixel) [0.05] - Size of the X axis in pixels [20] - Size of the Y axis in pixels [20] - Output Test Statistic map file [tsmap.fits] - - diff --git a/doc/source/users/tutorials/howto/iact/index.rst b/doc/source/users/tutorials/howto/iact/index.rst index a06cccf..006b36e 100755 --- a/doc/source/users/tutorials/howto/iact/index.rst +++ b/doc/source/users/tutorials/howto/iact/index.rst @@ -1,14 +1,17 @@ .. _sec_iact: -Analysing current IACT data -=========================== +Analysing data from existing IACTs +================================== + +This tutorial is specifically written for the analysis of H.E.S.S. data but +it may also be applicable to VERITAS and MAGIC, provided that these instruments +set up a data store that is compliant with the format described in +`Data formats for gamma-ray astronomy `_. .. toctree:: :maxdepth: 1 copy + select analysis - high_level - simulation - tips diff --git a/doc/source/users/tutorials/howto/iact/select.rst b/doc/source/users/tutorials/howto/iact/select.rst new file mode 100755 index 0000000..069db1e --- /dev/null +++ b/doc/source/users/tutorials/howto/iact/select.rst @@ -0,0 +1,224 @@ +.. _sec_iact_selection: + +Selecting IACT observations +=========================== + + .. admonition:: What you will learn + + You will learn how to select observations from a dataset. + +As next step you have to select observations from a dataset and to build +an observation definition XML file for the analysis. + + +Selecting observations +---------------------- + +First run :ref:`csfindobs` to select the observations of interest. +The most common search for observations is by pointing direction on the sky, +and you do this by typing + +.. code-block:: bash + + $ csfindobs + Path were data are located (NONE uses $VHEFITS environment variable) [NONE] data + Name of FITS production (Run csiactdata to view your options) [prod-name] hess_dl3_dr1 + Right ascension of selection region centre (deg) [83.63] + Declination of selection region centre (deg) [22.51] + Search radius of selection region (deg) [2.5] + Output runlist file [runlist.lis] + +This creates an output ASCII file named ``runlist.lis`` that contains the +identifiers for all observations with pointing directions located within 2.5 deg +around the position of the Crab nebula. The content of the file is + +.. code-block:: none + + 23523 + 23526 + 23559 + 23592 + +hence four observations were found in the production. An excerpt of the log +file ``csfindobs.log`` of the run is shown below: + +.. code-block:: none + + 2019-04-05T13:22:08: +===================+ + 2019-04-05T13:22:08: | Find observations | + 2019-04-05T13:22:08: +===================+ + 2019-04-05T13:22:08: Expression ................: ANGSEP(83.63,22.51,RA_PNT,DEC_PNT)<=2.5&&QUALITY<=0 + 2019-04-05T13:22:08: Observations ..............: 4 + 2019-04-05T13:22:08: Observation 1 .............: 23523 + 2019-04-05T13:22:08: Observation 2 .............: 23526 + 2019-04-05T13:22:08: Observation 3 .............: 23559 + 2019-04-05T13:22:08: Observation 4 .............: 23592 + +.. note:: + You may add further selection criteria using the hidden ``expression`` + parameter. + You may use any FITS selection expression that is + `supported by the cfitsio row filtering specification `_. + Column names that can be used in the expression can be found + `here `_ + or by simply browsing the observation index file with a FITS viewer. + For example + + .. code-block:: bash + + $ csfindobs expression="ZEN_PNT < 30 && LIVETIME > 1600" + Path were data are located (NONE uses $VHEFITS environment variable) [data] + Name of FITS production (Run csiactdata to view your options) [hess_dl3_dr1] + Right ascension of selection region centre (deg) [83.63] NONE + Output runlist file [runlist.lis] + + will select all observations that have a zenith angle less than 30 deg and a + livetime larger than 1600 seconds. Note that pointing selection was omitted + by specifying ``NONE`` for the Right Ascension. + +.. note:: + By default, :ref:`csfindobs` only selects data of highest quality + (i.e. ``QUALITY=0``). + You may overwrite this default by specifying the hidden parameter + ``min_qual``. For example, ``min_qual=1`` selects all data with a + looser quality criteria. + + +Create an observation list +-------------------------- + +Now you can convert the runlist into an observation definition XML file. +You do this conversion with the :ref:`csiactobs` script: + +.. code-block:: bash + + $ csiactobs + Path were data are located (NONE uses $VHEFITS environment variable) [NONE] data + Data storage name [fits-prod-name] hess_dl3_dr1 + Input runlist file [runlist.lis] + Number of free parameters per background model [1] 2 + Output model definition XML file [bkgmodels.xml] + Output observation definition XML file [obs.xml] + +The :ref:`csiactobs` script will create two output files: the observation +definition XML file ``obs.xml`` and an output model definition XML file +``bkgmodels.xml``. To generate ``obs.xml``, :ref:`csiactobs` has used the +IACT data storage and extracted the relevant file names. ``bgmodels.xml`` is +a file that is used for background modeling, where each observation will have +its own independent background model. In the example above, you have set the +number of free parameters per background model to one, hence the normalisation +of the background model for each observation will be a free parameter that is +later adjusted by the maximum likelihood fit. + +Here the observation definition file that was created + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +and here the model definition file + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +.. note:: + In case that you have already a model definition XML file that describes a + model of the celestial source distribution (a.k.a. a sky model), you may + provide this sky model to :ref:`csiactobs` using the hidden ``inmodel`` + parameter. The output model definition XML file will then contain the sky + model and the background model, and you can use the model directly in a + maximum likelihood fit. + + .. code-block:: bash + + $ csiactobs inmodel=$CTOOLS/share/models/crab.xml + Path were data are located (NONE uses $VHEFITS environment variable) [data] + Data storage name [hess_dl3_dr1] + Input runlist file [runlist.lis] + Number of free parameters per background model [2] + Output model definition XML file [bkgmodels.xml] models.xml + Output observation definition XML file [obs.xml] + + +Select events +------------- + +The final step needed is the selection of events, to define in particular +the region of interest and energy range for the analysis. +You do this by typing + +.. code-block:: bash + + $ ctselect + Input event list or observation definition XML file [events.fits] obs.xml + 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] + Lower energy limit (TeV) [0.1] 0.5 + Upper energy limit (TeV) [100.0] 10.0 + Output event list or observation definition XML file [selected_events.fits] obs_selected.xml + +which will select all events within a radius of 2 degrees around the pointing +direction with energies between 500 GeV and 10 TeV. + diff --git a/doc/source/users/tutorials/howto/iact/simulation.rst b/doc/source/users/tutorials/howto/iact/simulation.rst deleted file mode 100755 index 9e2ac45..0000000 --- a/doc/source/users/tutorials/howto/iact/simulation.rst +++ /dev/null @@ -1,300 +0,0 @@ -.. _sec_iact_simulation: - -Simulating IACT data -==================== -This section describes how to run simulations with current IACT data. In particular, a guide how to simulate data for a given set of IRFs, or a given set of -observations is provided. In addition, tools to compute sensitivity curves or pull distributions are presented. -All tools are based on the event list simulator :ref:`ctobssim`. - -Simulations for a specific observation container ------------------------------------------------- -On default, :ref:`ctobssim` builds a new observation container by querying relevant user parameters. In case -an observation XML container is already in place, the hidden parameter ``inobs`` can be specified (:ref:`see here how to create an -observation container `). It might be useful to run :ref:`ctselect` beforehand to ensure that the simulation knows about -energy thresholds, field of view, etc. If an unselected observation container is used, the simulation tool will query for global energy -boundaries and a global RoI radius. The following example illustrates the usage of an already selected observation container. - -.. code-block:: bash - - $ ctobssim inobs=selected_obs.xml - Input model XML file [crab_models.xml] - Output event data file or observation definition XML file [sim_obs.xml] - -Note that :ref:`ctobssim` writes a new event list file per observation into the current directory. If the observation container -is large, it might be convenient to create a new folder and specify the ``prefix`` keyword (similar to :ref:`ctselect`): - -.. code-block:: bash - - $ mkdir simulated - $ ctobssim inobs=selected_obs.xml prefix=simulated/sim_events_ - ... - -The newly created observation container can be used in the same way as a normal observation container, i.e. it can be analysed with -:ref:`ctlike` etc. - -.. note:: - - On default, :ref:`ctobssim` uses the same seed in every execution, i.e. the tool will always return the same simulated event list, regardless - when, where and how it was executed. To steer this the user can provide the hidden integer parameter ``seed``: - - .. code-block:: bash - - $ ctobssim inobs=selected_obs.xml seed=2 - - This might be important if simulations are run in parallel. - -Simulating observation time for a given set of IRFs ---------------------------------------------------- -There are other use cases to run simulations, for instance when there is no observation container available yet. This might be the case e.g -for an observation proposal: one would want to use a certain set of IRFs to simulate e.g. 50 hours of data to check whether a source -of interest can be detected. For this purpose it is necessary to first find a set of suitable IRFs to be used for -simulations. This can be easily done using :ref:`csfindobs` and using the ``expression`` parameter. The following example extracts -observations with at a certain zenith angle range (20-25 degrees) and azimuth angle range (90-100 degrees): - -.. code-block:: bash - - $ csfindobs expression="AZ_PNT>90.0&&AZ_PNT<100.0&&ZEN_PNT>20.0&&ZEN_PNT<25.0" - Name of FITS production (Run csiactdata to view your options) [fits-prod-name] - Right ascension [83.6331] NONE - Runlist outfile [zen20_az90.lis] - -Of course, further details could be specified, e.g. one might want to check for a certain range of muon efficieny, or any other parameter which impacts IRFs. -The output runlist subsequently has to be converted to an observation XML file using :ref:`csiactobs` (:ref:`see here how to execute this tool `). -In the following it is assumed that the observation XML file of the selected representative observations is produced and called ``myobs.xml``. - -Create caldb entry -^^^^^^^^^^^^^^^^^^ -Now, one of the observations in this container have to be added to the local calibration database. This database is used -by ctools to extract IRFs for simulations. - -.. code-block:: bash - - $ csobs2caldb - Input observation definition file [myobs.xml] - Response output name (e.g. Zenith50) [zen20_az90] - -Note that the response name has been specified to "zen20_az90". Any name can be chosen which represents the observation which was inserted into the database. -On default, :ref:`csobs2caldb` uses the first observation in the provided XML file. This can be modifyied by specifying the hidden integer parameter index, e.g.: - -.. code-block:: bash - - $ csobs2caldb index=1 - -This would use the second entry instead of the first from the observation container to extract the IRFs. - -.. note:: - - The ``index`` parameter starts counting from zero. - -Check caldb entries -^^^^^^^^^^^^^^^^^^^ -In order to check if the calibration database has successfully been added, the tool :ref:`cscaldb` can be used to inspect the calibration database: - -.. code-block:: bash - - $ cscaldb debug=yes - ... - 2016-02-25T15:57:11: +==============+ - 2016-02-25T15:57:11: | Mission: cta | - 2016-02-25T15:57:11: +==============+ - 2016-02-25T15:57:11: === Response functions in database "hess" === - 2016-02-25T15:57:11: zen20_az90 - 2016-02-25T15:57:11: - 2016-02-25T15:57:11: === Response functions in database "prod2" === - 2016-02-25T15:57:11: North_0.5h - 2016-02-25T15:57:11: North_50h - 2016-02-25T15:57:11: North_5h - 2016-02-25T15:57:11: South_0.5h - 2016-02-25T15:57:11: South_50h - 2016-02-25T15:57:11: South_5h - ... - -Note that the observation of the actual instrument ("hess" in this example) was added to the CTA mission. This is mainly for structural and simplicity reasons -and might be improved eventually. Nevertheless, this caldb entry can now be used to simulate any observation duration. - -Simulate -^^^^^^^^ - -.. code-block:: bash - - $ ctobssim - RA of pointing (degrees) (0-360) [83.63] - Dec of pointing (degrees) (-90-90) [22.01] - Radius of FOV (degrees) (0-180) [2.5] - Start time (MET in s) [0.0] - End time (MET in s) [18000.0] - Lower energy limit (TeV) [0.5] - Upper energy limit (TeV) [50] - Calibration database [hess] - Instrument response function [zen20_az90] - Input model XML file [input_model.xml] - Output event data file or observation definition XML file [simulated_events.fits] - -Here, a total duration of 5 hours at the crab position was simulated using the newly created caldb entry. Note that the choice of the energy range is up to the user. -One should bear in mind that simulations only make sense in energy ranges where the IRFs are well-defined. The simulation -assumes one observation of defined duration here. Accordingly, only one FITS event list is written as output. The input XML model -must now contain only one background model and (of course) all sky models that are needed. Here is an example of an input model -e.g. for the Crab Nebula: - -.. code-block:: xml - - - - - - - - - - - - - - - - - - - - - - - -Note that since the caldb entry (:ref:`csobs2caldb`) has been produced under the mission "CTA", the event list is now considered a CTA observation (regardless of the input instrument). -Therefore, it is important to have "CTA" specified in the instrument attribute of the background model XML component. Of course the addition of any sky model can be -done without a problem. - -The output of this kind of simulation is a single FITS event list file that can be used for further analysis. -For instance one could run :ref:`ctlike` in the following way: - -.. code-block:: bash - - $ ctlike - Input event list, counts cube or observation definition XML file [simulated_events.fits] - Calibration database [hess] - Instrument response function [zen20_az90] - Input model XML file [input_model.xml] - Output model XML file [simulated_fit_results.xml] - -Note that :ref:`ctlike` automatically queries for the IRF database in case no full observation container is provided. Any other tool taking -unbinned observations as input will behave the same way. Now all means to run simulations with current IACTs are at hand. -The following sections will describe tools for running higher level simulation analysis. - -Create pull distributions -------------------------- -Pull distributions inspect the capability of the science tools to find source parameters in the data that were injected into the simulations. -The tool :ref:`cspull` will run several times the sequence of simulating and fitting the data. After each sequence, the pull value of the -fitted parameters is written into an ASCII file. The pull value ``p`` of parameter ``X`` is defined as follows: - -.. math:: - p = \frac{X_{\rm fitted} - X_{\rm true}}{\Delta X_{\rm fitted}} - -Where :math:`\Delta X_{\rm fitted}` is the derived error of the fitted parameter :math:`X_{\rm fitted}`. After several trials, the pull value should follow -a normal distribution. Knowing this, the tool may help to distinguish if a source is detectable and can be normally reconstructed. -The tool can either work on an input observation container or on the caldb entry that was created above. The number of trials is queried by the tool, too. - -Example: observation container -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. code-block:: bash - - $ cspull inobs=selected_obs.xml - Number of energy bins (0=unbinned) [0] - Input model XML file [crab_results.xml] - Output pull distribution file [pull.dat] - Number of trials [100] - -Note that the tool asks for the energy binning. It will use an unbinned analysis if zero bins are provided. Otherwise -it will query for binning parameters and will convert the simulated observations into a binned/stacked observation by running -:ref:`ctbin`, :ref:`ctexpcube`, :ref:`ctpsfcube` and :ref:`ctbkgcube` in order before fitting the simulated data with :ref:`ctlike`. - -Example: caldb entry -^^^^^^^^^^^^^^^^^^^^ - -.. code-block:: bash - - $ cspull - RA of pointing (deg) (0-360) [83.6331] - Dec of pointing (deg) (-90-90) [22.0145] - Duration (in s) [18000.0] - Lower energy limit (TeV) [0.5] - Upper energy limit (TeV) [50] - Calibration database [hess] - Instrument response function [zen20_az90] - Number of energy bins (0=unbinned) [0] - Input model XML file [input_models.xml] - Output pull distribution file [pull.dat] - Number of trials [100] - -Again, the caldb entry ``zen20_az90`` which was created above is used to simulate 5 hours of data. The data is simulated and analysed 100 -times in a row and the output is written into the file ``pull.dat``. - -Visualise pull distributions ----------------------------- -The output ASCII file of :ref:`cspull` can be inspected by the user using own scripts. As a starting point, the ``example`` folder -contains a script that can be used to visualise the pull distribution histogram: - -.. code-block:: bash - - $ python $CTOOLS/examples/show_pull_histogram.py pull.dat Pull_Crab_Prefactor 50 - -Run the script without arguments to see the usage. In case a wrong parameter name was provided, the tool will print available parameters on the screen. -The last integer parameter defines the binning of the histogram. In addition, the tool shows the standard normal distribution which should -match the pull distribution for a high number of trials if everything went well. - -Create sensitivity curve ------------------------- -A very important means to illustrate if a source can be detected with a certain significance is a sensitivity curve. -The tool :ref:`cssens` can be used to compute either an integral, or a differential sensitivity curve for a given observation container. -Analogous to :ref:`cspull`, the tool works either on an unbinned observation container or using a caldb entry. -The script :ref:`cssens` computes the differential or integrated CTA sensitivity using maximum likelihood fitting of a test source. -The differential sensitivity is determined for a number of energy bins, the integral sensitivity is determined for a number -of energy thresholds. The test source is fitted to simulated data using ctlike to determine it’s detection significance as a -function of source flux. The source flux is then varied until the source significance achieves a given level, specified by the -(hidden) significance parameter ``sigma``. To damp variations between individual Monte Carlo simulations, a sliding average is -applied in the significance computation (controlled by hidden ``num_avg parameter``). Note that this procedure might take a while to compute. -To switch from diffferential to integral sensitivity, the hidden parameter ``type="Integral"`` has to be specified. It should be noted that -the results are written to a file ``sensitivity.dat``. To change the output file name, the hidden parameter ``outfile=mysensitivity.dat`` -can be provided. - - -Example: observation container -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. code-block:: bash - - $ cssens inobs=selected_obs.xml - Input model XML file [crab_results.xml] - Source name [Crab] - Lower energy limit (TeV) [0.5] - Upper energy limit (TeV) [50] - Number of energy bins for differential sensitivity computation [21] - -Example: caldb entry -^^^^^^^^^^^^^^^^^^^^ - -.. code-block:: bash - - $ cssens - Calibration database [hess] - Instrument response function [zen20_az90] - Effective exposure time (s) [18000.0] - Radius of ROI (deg) [2.5] - Input model XML file [input_models.xml] - Source name [Crab] - Lower energy limit (TeV) [0.5] - Upper energy limit (TeV) [50] - Number of energy bins for differential sensitivity computation [21] - -Note that the input model in the latter case only requires to contain one background model (see section above). Using this setup, -the user can specify the observation time for the sensitivity curve. In this example, again, 5 hours are used. - -Visualise sensitivity curves ----------------------------- -In order to show a sensitivity curve, there is also a script in the example folder: - -.. code-block:: bash - - $ python $CTOOLS/examples/show_sensitivity.py sensitivity.dat - -This will display a sensitivity curve for the applied settings. diff --git a/doc/source/users/tutorials/howto/iact/tips.rst b/doc/source/users/tutorials/howto/iact/tips.rst deleted file mode 100755 index 1f09b91..0000000 --- a/doc/source/users/tutorials/howto/iact/tips.rst +++ /dev/null @@ -1,503 +0,0 @@ -.. _sec_iact_tips: - -Tips, Tricks and FAQs -===================== -This section shows some of the wrinkles that could ease your life when -running IACT analyses. - - -Visualise observations ----------------------- -A very common task is to inspect the observation container and visualise several -quantities. - -Create a skymap of counts -^^^^^^^^^^^^^^^^^^^^^^^^^ -The tool :ref:`ctskymap` allows to create a single map of counts contained -inside an observation container. In contrast to :ref:`ctbin`, ctskymap takes -into account all events that fall into the skymap energy range and spatial -range. Since :ref:`ctbin` is intended for stacked analysis, it respects energy -thresholds, RoI selections, etc. In contrast, :ref:`ctskymap` may be used for -a simple visualisation of the counts. - -.. code-block:: bash - - $ ctskymap debug=yes - Input event list or observation definition XML file [events.fits] obs.xml - Lower energy limit (TeV) [0.1] - Upper energy limit (TeV) [100.0] - 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] - Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] - Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] - Image scale (in degrees/pixel) [0.02] - Size of the X axis in pixels [200] - Size of the Y axis in pixels [200] - Output skymap file [skymap.fits] - -Alternatively, you may start from a counts cube and collapse the maps into one -single map via Python: - -.. code-block:: python - - >>> import gammalib - >>> map = gammalib.GSkyMap("myCountCube.fits") # Load counts cube - >>> map.stack_maps() # Stack all cube maps into a single map - >>> map.save("myCountMap.fits") # Save map to file - -Create region file of pointings -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -You may want to know how the pointings of the observation container distribute -across the sky. For this purpose, :ref:`csobsinfo` provides the hidden parameter -``ds9file`` which instructs the script produce a region file that contains the -pointing information: - -.. code-block:: bash - - $ csobsinfo ds9file=pointings.reg debug=yes - Event list, counts cube, or observation definition file [obs.xml] - -Together with the previously produced skymap, you can plot the pointings using -for example ds9: - -.. code-block:: bash - - $ ds9 skymap.fits -regions pointings.reg - - -Plot zenith angle distribution of observations ----------------------------------------------- -The script :ref:`csobsinfo` logs information about an observation container -into a log file. You may access the information through dedicated class methods -and plot the information using for example matplotlib: - -.. code-block:: python - - >>> # Run script - >>> import cscripts - >>> info = cscripts.csobsinfo() - >>> info["inobs"] = "selected_obs.xml" - >>> info.run() - - >>> # Plot data - >>> import matplotlib.pyplot as plt - >>> plt.hist(info.zeniths(), bins=10, range=[0,90]) - >>> plt.show() - -This example will plot a histogram showing the zenith angle distribution of the -observation container. In addition there are more methods to access information -in Python: - -.. code-block:: python - - >>> info.azimuths() # list of azimuth values - >>> info.offsets() # list of offset values (only computed if hidden parameter offset=yes was specified) - >>> info.ebounds() # gammalib.GEbounds object of energy ranges - >>> info.gti() # gammalib.GGti object containing good time intervals - -There is also an example script to plot information about an observation -container which might be useful for a quick visualisation: - -.. code-block:: bash - - $ python $CTOOLS/examples/show_obs.py selected_obs.xml - - -Visualise models ----------------- -Similar to :ref:`csobsinfo`, there is the :ref:`csmodelinfo` script to find out -more about the content of a model XML file. In order to show the position and -sizes of the model on top of a skymap, this tool has the hidden parameter -``ds9file``, too. - -.. code-block:: bash - - $ csmodelinfo ds9file=models.reg debug=yes - Input model XML file [$CTOOLS/share/models/crab.xml] - -There are several options regarding the color, text and other attributes of the -region file. To see a full list, :ref:`visit the reference page `. -Analogous to the pointings of the observation container, the models can be -visualised using for example ds9: - -.. code-block:: bash - - $ ds9 skymap.fits -regions models.reg - - -Manipulating models in python ------------------------------ -This example will give you some guidance on how to work with model XML files in -Python and how to manipulate their content: - -.. code-block:: python - - >>> import gammalib - - >>> # Open model file - >>> models = gammalib.GModels("$CTOOLS/share/models/crab.xml") - - >>> # Access a model component - >>> src = models["Crab"] - >>> print(src) - - >>> # Retrieve the spectral or spatial component - >>> print(src.spectral()) - >>> print(src.spatial()) - - >>> # Print the spectral parameter "Prefactor" - >>> prefactor = src.spectral()["Prefactor"] - >>> print("Prefactor: "+str(prefactor.value())+" +- "+str(prefactor.error())) - - >>> # Set the prefactor value - >>> prefactor.value(3.5e-17) - - >>> # Loop over models and fix all parameters - >>> for model in models: - >>> for par in model: # Loop over model parameters - >>> par.fix() # Fix parameter - - >>> # Release one specific parameter - >>> models["Crab"]["Prefactor"].free() - - >>> # Set the parameter range - >>> models["Crab"]["Prefactor"].min(1e-18) - >>> models["Crab"]["Prefactor"].max(1e-16) - - >>> # ... or in one step: - >>> models["Crab"]["Prefactor"].range(1e-18, 1e-16) - - >>> # Remove model from container - >>> models.remove("Crab") - - >>> # Save model to another XML file - >>> models.save("mymodels.xml") - - -Retrieve likelihood values from ctlike --------------------------------------- -For some purposes, it might be useful to retrieve the fitted likelihood value -and other results of the fit with :ref:`ctlike`. - -.. code-block:: python - - >>> import gammalib - >>> import ctools - - >>> # Create and run ctlike - >>> like = ctools.ctlike() - >>> like["inobs"] = "selected_obs.xml" - >>> like["inmodel"] = "mymodels.xml" - >>> like.run() - - >>> # Get optimizer from ctlike - >>> opt = like.opt() - >>> print(opt) - - >>> # Get likelihood value from optimizer - >>> print(opt.value()) - - >>> # Get fit status: - >>> print(opt.status()) - >>> # 0: converged - >>> # 1: stalled - >>> # 2: singular curvature matrix encountered - >>> # 3: curvature matrix not positive definite - >>> # 4: errors are inaccurate - - -Speed up analysis ------------------ -In some cases it may occur that the analysis takes very long. There are several -reasons that can slow down the analysis quite severly. Not all issues can be -tackled. Nevertheless, here is a list of actions that have proven to speed up -the fit. - -* Consider switching to stacked analysis if the observation time exceeds 50-100 - hours. -* Diffuse and extended sources in the RoI are computational expensive. However, - in unbinned analysis they impact the speed much more. -* Usage of energy dispersion could also cause a long fit. Try switching it off - if speed is important. Keep in mind this will have an effect on the results. -* Switch to a machine that supports OpenMP. This allows to parallelise the fit - onto several cores. - - -Compute excess maps -------------------- -A very important means to visualise images of sources are excess maps. The script -:ref:`csresmap` can also be used to create such a map. The user, however, has to -manually remove the source(s) of interest from the input model XML file. - -.. code-block:: python - - >>> import gammalib - >>> models = gammalib.GModels("crab_models.xml") - >>> models.remove("Crab") - >>> models.save("models_without_crab.xml") - -Of course, this can also be done by editing the XML file with an editor. - -Subsequently, :ref:`csresmap` can be executed using the ``algorithm=SUB`` -parameter: - -.. code-block:: bash - - $ csresmap - Input event list, counts cube, or observation definition XML file [selected_obs.xml] - Input model XML file [models_without_crab.xml] - 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] - Coordinate System (CEL|GAL) [CEL] - Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] - Size of the X axis in pixels [100] - Size of the Y axis in pixels [100] - Pixel size (deg/pixel) [0.02] - Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUB] - Output residual map file [excessmap.fits] - - -Compute flux maps ------------------ -The creation of a flux map is much more expensive in terms of computing than an -excess map. In the section about :ref:`high level analysis tools `, -the usage of :ref:`cttsmap` was demonstrated. This tool can also be used to -generate flux maps. For this, the user has to incorporate a test source in the -XML model file. This source should consist of a spectral component using a -``PowerLawPhotonFlux`` model: - -.. code-block:: xml - - - - - - - - - - - - - - -The further content of the rest of the XML model file depends on the user -requirements: - -* For an absolute flux map, all other sky models should be removed -* For a residual flux map, all other sky models should be kept - -The tool :ref:`cttsmap` will create one skymap per free spectral parameter in -the model. In the result FITS file, there will be an extension called ``Integral`` -that contains the flux map. For the above XML example, the source name -``TestSource`` should be specified to :ref:`cttsmap`. - - -Speed up TS map computation ---------------------------- -Since the computation of a TS map can be extremely time consuming, the option of -splitting the computation into several jobs is supported. This might be of -particular interest if the user has access to a batch farm. - -Split TS map computation -^^^^^^^^^^^^^^^^^^^^^^^^ -For the purpose of job splitting, the hidden parameter ``binmin`` and ``binmax`` -were included in the tool. These are integer parameter that specify which bins -should be computed. For instance, if the map should consist of 30x30(=900) pixels, -the user could for example execute - -.. code-block:: bash - - $ cttsmap binmin=0 binmax=299 outmap=tsmap_0_299.fits - ... - $ cttsmap binmin=300 binmax=599 outmap=tsmap_300_599.fits - ... - $ cttsmap binmin=600 binmax=899 outmap=tsmap_600_899.fits - ... - -Each command could run on a different machine. - -.. note:: - - The output file name of individual jobs should be different. Otherwise files - could overwrite each other. The naming of the individual slices is up to the - user. - -The script :ref:`cstsmapsplit` will take care of this bookkeeping. It creates an -ASCII file containing all the commands according to the user input. The following -example will create split the computation of a TS map into 2000 separate task. -Each task will only compute 5 bins on its own. This is very useful if the -observation container is large or the fit simply takes a long time. - -.. code-block:: bash - - $ cstsmapsplit - Input event list, counts cube or observation definition XML file [selected_obs.xml] - Input model definition XML file [$CTOOLS/share/models/crab.xml] - 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] - Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] - Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] - Image scale (in degrees/pixel) [0.02] - Size of the X axis in pixels [100] - Size of the Y axis in pixels [100] - Test source name [Crab] - Output Test Statistic map file [tsmap.fits] - Number of TS map bins per task [5] - Compute null hypothesis first? [yes] - ASCII file containing all commands [commands.dat] - -.. note:: - - One can decide if we want to compute the null hypothesis first. This way, a - parameter optimisation will be performed without the test source. The obtained - likelihood value will then be passed to each individual task of :ref:`cttsmap`. - If ``compute_null=no``, each task has to compute the null hypothesis itself. - -Merge splitted TS maps -^^^^^^^^^^^^^^^^^^^^^^ -The script :ref:`cstsmapmerge` is intended to take care of merging the files -that were produced while splitting the TS map computation. There are several -options to pass the files to be merged as arguments: - -* a space-separated list of file names (e.g. tsmap1.fits tsmap2.fits) -* a comma-separated list of file names (e.g. tsmap1.fits,tsmap2.fits) -* a wildcard string (e.g. tsmap*.xml) -* an ASCII file containing the file names, one file per line (e.g. @mymaps.txt) - -In this example the ASCII file method is presented: - -.. code-block:: bash - - # Put slice files into an ascii file - $ ls tsmap_*.fits > tsmapfiles.txt - - # Run cstsmapmerge - $ cstsmapmerge - Input TS map FITS files [@tsmapfiles.txt] - Output TS map FITS file [mytsmap.fits] - - -Creating a python analysis pipeline ------------------------------------ -It is easily possible to build an own analysis workflow with a simple python -script. The following source code example shows a python script running from -gathering observations until fitting spectral points without storing intermediate -data products on disk. It assumes that the environment variable ``$VHEFITS`` is -set to the path where IACT FITS data is located. - -.. code-block:: python - - >>> import gammalib - >>> import ctools - >>> import cscripts - - >>> # Set debug flag - >>> debug = True - - >>> # Set flag to use energy dispersion - >>> edisp = False - - >>> # Set inmodel file name - >>> inmodel = "$GAMMALIB/test/data/model_point_plaw.xml" - - >>> # Expand environment variable - >>> inmodel = gammalib.expand_env(inmodel) - - >>> # Extract coordinates and model properties - >>> models = gammalib.GModels(inmodel) - >>> srcname = models[0].name() - >>> ra = models[0]["RA"].value() - >>> dec = models[0]["DEC"].value() - - >>> # Find FITS production name - >>> iactdata = cscripts.csiactdata() - >>> iactdata["debug"] = debug - >>> iactdata.run() - - >>> # Use first available production - >>> prodname = iactdata.names()[0] - - >>> # Run csfindobs - >>> findobs = cscripts.csfindobs() - >>> findobs["ra"] = ra - >>> findobs["dec"] = dec - >>> findobs["rad"] = 2.5 - >>> findobs["prodname"] = prodname - >>> findobs["debug"] = debug - >>> findobs["outfile"] = "NONE" - >>> findobs.run() - - >>> # Retrieve obervation IDs (runlist) - >>> obs_ids = findobs.obs_ids() - - >>> # Build observation container - >>> iactobs = cscripts.csiactobs() - >>> iactobs["prodname"] = prodname - >>> iactobs["inmodel"] = inmodel - >>> iactobs["bkgpars"] = 1 - >>> iactobs["outobs"] = "NONE" - >>> iactobs["outmodel"] = "NONE" - >>> iactobs["debug"] = debug - >>> iactobs.runlist(obs_ids) - >>> iactobs.run() - - >>> # Retrieve observation container and energy boundaries - >>> obs = iactobs.obs() - >>> ebounds = iactobs.ebounds() - - >>> # Run ctselect - >>> select = ctools.ctselect(obs) - >>> select["usepnt"] = True - >>> select["rad"] = 2.5 - >>> select["usethres"] = "DEFAULT" - >>> select["tmin"] = 0.0 - >>> select["tmax"] = 0.0 - >>> select["emin"] = 0.1 - >>> select["emax"] = 100.0 - >>> select["debug"] = debug - >>> select.run() - - >>> # Pass selected observations to ctlike - >>> like = ctools.ctlike(select.obs()) - >>> like["debug"] = debug - >>> like["edisp"] = edisp - >>> like.run() - - >>> # Compute a spectrum and save - >>> spec = cscripts.csspec(like.obs()) - >>> spec["srcname"] = srcname - >>> spec["emin"] = ebounds.emin().TeV() - >>> spec["emax"] = ebounds.emax().TeV() - >>> spec["enumbins"] = 10 - >>> spec["edisp"] = edisp - >>> spec["ebinalg"] = "LOG" - >>> spec["debug"] = debug - >>> spec["outfile"] = "spectrum.fits" - >>> spec.execute() - - >>> # Remove model for excess map computation - >>> like.obs().models().remove(srcname) - - >>> # Compute an excess map - >>> resmap = cscripts.csresmap(like.obs()) - >>> resmap["xref"] = ra - >>> resmap["yref"] = dec - >>> resmap["proj"] = "CAR" - >>> resmap["coordsys"] = "CEL" - >>> resmap["emin"] = ebounds.emin().TeV() - >>> resmap["emax"] = ebounds.emax().TeV() - >>> resmap["nxpix"] = 100 - >>> resmap["nypix"] = 100 - >>> resmap["binsz"] = 0.02 - >>> resmap["outmap"] = "excessmap.fits" - >>> resmap.execute() - -The results of this analysis workflow can be inspected using an example Python -script and for example ds9: - -.. code-block:: bash - - $ python $CTOOLS/examples/show_spectrum.py spectrum.fits - $ ds9 excessmap.fits - diff --git a/doc/source/users/tutorials/howto/index.rst b/doc/source/users/tutorials/howto/index.rst index a8b22ab..bba0870 100755 --- a/doc/source/users/tutorials/howto/index.rst +++ b/doc/source/users/tutorials/howto/index.rst @@ -46,5 +46,5 @@ Analysis of data from other instruments How to analyse Fermi-LAT data with ctools? How to analyse COMPTEL data with ctools? - How to analyse H.E.S.S. data with ctools? + How to analyse data from existing IACTs (H.E.S.S., VERITAS or MAGIC) with ctools? howto_combine_instruments