diff --git a/doc/source/users/tutorials/howto/fermi/howto_fermi_butterfly.png b/doc/source/users/tutorials/howto/fermi/howto_fermi_butterfly.png new file mode 100644 index 0000000..007abed --- /dev/null +++ b/doc/source/users/tutorials/howto/fermi/howto_fermi_butterfly.png diff --git a/doc/source/users/tutorials/howto/fermi/howto_fermi_butterfly.rst b/doc/source/users/tutorials/howto/fermi/howto_fermi_butterfly.rst new file mode 100644 index 0000000..ca26b05 --- /dev/null +++ b/doc/source/users/tutorials/howto/fermi/howto_fermi_butterfly.rst @@ -0,0 +1,48 @@ +.. _howto_fermi_butterfly: + +Generate a butterfly diagram for a model component +-------------------------------------------------- + + .. admonition:: What you will learn + + You will learn how to use :ref:`ctbutterfly` to create a butterfly diagram + for a model component from Fermi-LAT data. + +A butterfly diagram shows the envelope of all spectral models that are +statistically compatible with the data. You can create such a butterfly +diagram for the Vela pulsar using the :ref:`ctbutterfly` tool by providing +the +:ref:`model definition file ` +output file from :ref:`ctlike` as input: + +.. code-block:: bash + + $ ctbutterfly + Input event list, counts cube or observation definition XML file [events.fits] obs.xml + Source of interest [Crab] Vela + Input model definition XML file [$CTOOLS/share/models/crab.xml] vela_results.xml + Start value for first energy bin in TeV [0.1] 0.0001 + Stop value for last energy bin in TeV [100.0] 0.1 + Output ASCII file [butterfly.txt] + +:ref:`ctbutterfly` produces the ASCII file ``butterfly.txt`` that contains +the envelope of the butterfly diagram. The figure below gives a graphical +representation of the envelope. + +.. figure:: howto_fermi_butterfly.png + :width: 600px + :align: center + + *Butterfly diagram of the fitted Vela pulsar spectrum* +.. + + .. note:: + + The figure was created using the ``show_butterfly.py`` script that is + located in the ctools example folder. The example script requires the + `matplotlib `_ Python module for display. + You may reproduce the plot by typing + + .. code-block:: bash + + $ $CTOOLS/share/examples/python/show_butterfly.py butterfly.txt diff --git a/doc/source/users/tutorials/howto/fermi/howto_fermi_data.jpg b/doc/source/users/tutorials/howto/fermi/howto_fermi_data.jpg new file mode 100644 index 0000000..bd17969 --- /dev/null +++ b/doc/source/users/tutorials/howto/fermi/howto_fermi_data.jpg diff --git a/doc/source/users/tutorials/howto/fermi/howto_fermi_errors.rst b/doc/source/users/tutorials/howto/fermi/howto_fermi_errors.rst new file mode 100644 index 0000000..79dadd4 --- /dev/null +++ b/doc/source/users/tutorials/howto/fermi/howto_fermi_errors.rst @@ -0,0 +1,86 @@ +.. _howto_fermi_errors: + +Determine asymmetric model parameter errors for a model component +----------------------------------------------------------------- + + .. admonition:: What you will learn + + You will learn how to use :ref:`cterror` to derive asymmetric model + parameter errors from Fermi-LAT data using a likelihood profile method. + +The :ref:`ctlike` tool computes errors for each fitted parameter value based +on the covariance matrix of the fit. This assumes a Gaussian approximation +of the uncertainties around the best fitting values. Uncertainties may however +by asymmetric. + +To test the symmetry of the uncertainties you can use the :ref:`cterror` +tool. :ref:`cterror` will make use of the likelihood profile to determine +the confidence interval for each fitted parameter, allowing thus also for +asymmetric errors. You use :ref:`cterror` as follows: + +.. code-block:: bash + + $ cterror + Input event list, counts cube or observation definition XML file [events.fits] obs.xml + Source of interest [Crab] Vela + Input model definiton XML file [$CTOOLS/share/models/crab.xml] vela_results.xml + Output model definiton XML file [results.xml] vela_errors.xml + +The confidence intervals can then be extracted from the ``cterror.log`` log file +that was created by :ref:`cterror`. As you will see, the parameter errors +are essentially symmetric and very close to the values determined using the +covariance matrix. + +.. code-block:: bash + + 2017-08-26T14:16:00: +=======================================================+ + 2017-08-26T14:16:00: | Compute error for source "Vela" parameter "Prefactor" | + 2017-08-26T14:16:00: +=======================================================+ + 2017-08-26T14:16:00: Confidence level ..........: 68 % + 2017-08-26T14:16:00: Log-likelihood difference .: 0.494473 + 2017-08-26T14:16:00: Initial factor range ......: [1e-07, 9.10916268820218] + 2017-08-26T14:16:20: Lower parameter factor ....: 3.93549 + 2017-08-26T14:16:20: Upper parameter factor ....: 4.89285 + 2017-08-26T14:16:20: Error from curvature ......: 4.75365104151534e-10 ph/cm2/s/MeV + 2017-08-26T14:16:20: Error from profile ........: 4.78682555317737e-10 ph/cm2/s/MeV + 2017-08-26T14:16:20: Negative profile error ....: 4.20026137925122e-10 ph/cm2/s/MeV + 2017-08-26T14:16:20: Positive profile error ....: 5.37338972710351e-10 ph/cm2/s/MeV + 2017-08-26T14:16:20: + 2017-08-26T14:16:20: +====================================================+ + 2017-08-26T14:16:20: | Compute error for source "Vela" parameter "Index1" | + 2017-08-26T14:16:20: +====================================================+ + 2017-08-26T14:16:20: Confidence level ..........: 68 % + 2017-08-26T14:16:20: Log-likelihood difference .: 0.494473 + 2017-08-26T14:16:20: Initial factor range ......: [1.03221907174349, 1.65479046234989] + 2017-08-26T14:16:42: Lower parameter factor ....: 1.31121 + 2017-08-26T14:16:42: Upper parameter factor ....: 1.37307 + 2017-08-26T14:16:42: Error from curvature ......: 0.03112856953032 + 2017-08-26T14:16:42: Error from profile ........: 0.0309309760713872 + 2017-08-26T14:16:42: Negative profile error ....: 0.0322989307870749 + 2017-08-26T14:16:42: Positive profile error ....: 0.0295630213556994 + 2017-08-26T14:16:42: + 2017-08-26T14:16:42: +==========================================================+ + 2017-08-26T14:16:42: | Compute error for source "Vela" parameter "CutoffEnergy" | + 2017-08-26T14:16:42: +==========================================================+ + 2017-08-26T14:16:42: Confidence level ..........: 68 % + 2017-08-26T14:16:42: Log-likelihood difference .: 0.494473 + 2017-08-26T14:16:42: Initial factor range ......: [0.001, 2.8071605417966] + 2017-08-26T14:17:02: Lower parameter factor ....: 0.813181 + 2017-08-26T14:17:02: Upper parameter factor ....: 1.17404 + 2017-08-26T14:17:02: Error from curvature ......: 181.871842544747 MeV + 2017-08-26T14:17:02: Error from profile ........: 180.431484426618 MeV + 2017-08-26T14:17:02: Negative profile error ....: 175.261332662553 MeV + 2017-08-26T14:17:02: Positive profile error ....: 185.601636190684 MeV + 2017-08-26T14:17:02: + 2017-08-26T14:17:02: +====================================================+ + 2017-08-26T14:17:02: | Compute error for source "Vela" parameter "Index2" | + 2017-08-26T14:17:02: +====================================================+ + 2017-08-26T14:17:02: Confidence level ..........: 68 % + 2017-08-26T14:17:02: Log-likelihood difference .: 0.494473 + 2017-08-26T14:17:02: Initial factor range ......: [0.307804018748382, 0.866425033775455] + 2017-08-26T14:17:24: Lower parameter factor ....: 0.559156 + 2017-08-26T14:17:24: Upper parameter factor ....: 0.615141 + 2017-08-26T14:17:24: Error from curvature ......: 0.0279310507513537 + 2017-08-26T14:17:24: Error from profile ........: 0.0279924226890397 + 2017-08-26T14:17:24: Negative profile error ....: 0.027958327168103 + 2017-08-26T14:17:24: Positive profile error ....: 0.0280265182099764 diff --git a/doc/source/users/tutorials/howto/fermi/howto_fermi_fitting.rst b/doc/source/users/tutorials/howto/fermi/howto_fermi_fitting.rst new file mode 100644 index 0000000..833b206 --- /dev/null +++ b/doc/source/users/tutorials/howto/fermi/howto_fermi_fitting.rst @@ -0,0 +1,256 @@ +.. _howto_fermi_fitting: + +Fit a model to the data +----------------------- + + .. admonition:: What you will learn + + You will learn how to use :ref:`ctlike` to fit a parametric model to + Fermi-LAT data. + +Start with creating an +:ref:`observation definition file ` +that defines the names of the counts map file, the exposure map file, the +live time cube file and the name of the instrument response function that +should be used for the analysis: + +.. code-block:: xml + + + + + + + + + + + +This will be your input file to all tools. +Then you also need a +:ref:`model definition file ` +that describes the source and background model that you want to fit. In the +example below, a point source with a super exponentially cut-off power law is +fit at the position of the Vela pulsar on top of a Galactic diffuse and an +extragalactic isotropic source component. + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +.. + + .. warning:: + + The names of the diffuse model components have to be identical to the + names used when preparing the ``srcmaps.fits`` file using ``gtsrcmap``. + + **Positions of point sources can so far not be fitted for a Fermi-LAT + analysis**, hence the corresponding parameters are fixed. + +Now you are ready to fit the model using the :ref:`ctlike` model tool. You do +this as follows: + +.. code-block:: bash + + $ ctlike + Input event list, counts cube or observation definition XML file [events.fits] obs.xml + Input model definition XML file [$CTOOLS/share/models/crab.xml] models.xml + Output model definition XML file [crab_results.xml] vela_results.xml + +:ref:`ctlike` creates on output the +:ref:`model definition file ` +``vela_results.xml`` where the parameter values were updated by their +fitted values, and where the statistical uncertainties were added using the +``error`` attribute. To investigate how the fit went you should inspect the +log file ``ctlike.log`` that was also created by :ref:`ctlike`: + +.. code-block:: bash + + 2017-08-26T10:35:10: +=================================+ + 2017-08-26T10:35:10: | Maximum likelihood optimisation | + 2017-08-26T10:35:10: +=================================+ + 2017-08-26T10:35:12: >Iteration 0: -logL=-177232.327, Lambda=1.0e-03 + 2017-08-26T10:35:12: >Iteration 1: -logL=-198931.081, Lambda=1.0e-03, delta=21698.754, step=1.0e+00, max(|grad|)=-22713.401202 [CutoffEnergy:4] + 2017-08-26T10:35:12: >Iteration 2: -logL=-200263.519, Lambda=1.0e-04, delta=1332.438, step=1.0e+00, max(|grad|)=-67648.133552 [CutoffEnergy:4] + 2017-08-26T10:35:12: >Iteration 3: -logL=-212350.035, Lambda=1.0e-05, delta=12086.515, step=1.0e+00, max(|grad|)=23036.151920 [Index2:6] + 2017-08-26T10:35:12: >Iteration 4: -logL=-213951.781, Lambda=1.0e-06, delta=1601.746, step=1.0e+00, max(|grad|)=-8297.648633 [CutoffEnergy:4] + 2017-08-26T10:35:13: >Iteration 5: -logL=-215341.362, Lambda=1.0e-07, delta=1389.581, step=1.0e+00, max(|grad|)=2747.034517 [Index2:6] + 2017-08-26T10:35:13: >Iteration 6: -logL=-215490.429, Lambda=1.0e-08, delta=149.067, step=1.0e+00, max(|grad|)=697.301361 [Index2:6] + 2017-08-26T10:35:13: Iteration 7: -logL=-215490.429, Lambda=1.0e-09, delta=-14.065, step=1.0e+00, max(|grad|)=-1205.683279 [Index1:3] (stalled) + 2017-08-26T10:35:13: Iteration 8: -logL=-215490.429, Lambda=1.0e-08, delta=-14.064, step=1.0e+00, max(|grad|)=-1205.659884 [Index1:3] (stalled) + 2017-08-26T10:35:14: Iteration 9: -logL=-215490.429, Lambda=1.0e-07, delta=-14.055, step=1.0e+00, max(|grad|)=-1205.425967 [Index1:3] (stalled) + 2017-08-26T10:35:14: Iteration 10: -logL=-215490.429, Lambda=1.0e-06, delta=-13.963, step=1.0e+00, max(|grad|)=-1203.090620 [Index1:3] (stalled) + 2017-08-26T10:35:14: Iteration 11: -logL=-215490.429, Lambda=1.0e-05, delta=-13.066, step=1.0e+00, max(|grad|)=-1180.113526 [Index1:3] (stalled) + 2017-08-26T10:35:14: Iteration 12: -logL=-215490.429, Lambda=1.0e-04, delta=-6.263, step=1.0e+00, max(|grad|)=-983.079153 [Index1:3] (stalled) + 2017-08-26T10:35:14: >Iteration 13: -logL=-215496.367, Lambda=1.0e-03, delta=5.937, step=1.0e+00, max(|grad|)=-285.734284 [Index1:3] + 2017-08-26T10:35:15: Iteration 14: -logL=-215496.367, Lambda=1.0e-04, delta=-8.812, step=1.0e+00, max(|grad|)=-731.894749 [Index1:3] (stalled) + 2017-08-26T10:35:15: >Iteration 15: -logL=-215497.536, Lambda=1.0e-03, delta=1.169, step=1.0e+00, max(|grad|)=-153.668138 [Index1:3] + 2017-08-26T10:35:15: Iteration 16: -logL=-215496.136, Lambda=1.0e-04, delta=-1.400, step=1.0e+00, max(|grad|)=-328.522942 [Index1:3] (stalled) + 2017-08-26T10:35:15: >Iteration 17: -logL=-215498.376, Lambda=1.0e-03, delta=2.240, step=1.0e+00, max(|grad|)=3.257091 [Index2:6] + 2017-08-26T10:35:16: >Iteration 18: -logL=-215498.383, Lambda=1.0e-04, delta=0.007, step=1.0e+00, max(|grad|)=-8.270351 [CutoffEnergy:4] + 2017-08-26T10:35:16: >Iteration 19: -logL=-215498.384, Lambda=1.0e-05, delta=0.002, step=1.0e+00, max(|grad|)=0.637485 [Index2:6] + 2017-08-26T10:35:16: + 2017-08-26T10:35:16: +=========================================+ + 2017-08-26T10:35:16: | Maximum likelihood optimisation results | + 2017-08-26T10:35:16: +=========================================+ + 2017-08-26T10:35:16: === GOptimizerLM === + 2017-08-26T10:35:16: Optimized function value ..: -215498.384 + 2017-08-26T10:35:16: Absolute precision ........: 0.005 + 2017-08-26T10:35:16: Acceptable value decrease .: 2 + 2017-08-26T10:35:16: Optimization status .......: converged + 2017-08-26T10:35:16: Number of parameters ......: 14 + 2017-08-26T10:35:16: Number of free parameters .: 5 + 2017-08-26T10:35:16: Number of iterations ......: 19 + 2017-08-26T10:35:16: Lambda ....................: 1e-06 + 2017-08-26T10:35:16: Maximum log likelihood ....: 215498.384 + 2017-08-26T10:35:16: Observed events (Nobs) ...: 202330.000 + 2017-08-26T10:35:16: Predicted events (Npred) ..: 202318.082 (Nobs - Npred = 11.9178772509331) + 2017-08-26T10:35:16: === GModels === + 2017-08-26T10:35:16: Number of models ..........: 3 + 2017-08-26T10:35:16: Number of parameters ......: 14 + 2017-08-26T10:35:16: === GModelSky === + 2017-08-26T10:35:16: Name ......................: Vela + 2017-08-26T10:35:16: Instruments ...............: all + 2017-08-26T10:35:16: Instrument scale factors ..: unity + 2017-08-26T10:35:16: Observation identifiers ...: all + 2017-08-26T10:35:16: Model type ................: PointSource + 2017-08-26T10:35:16: Model components ..........: "PointSource" * "SuperExponentialCutoffPowerLaw" * "Constant" + 2017-08-26T10:35:16: Number of parameters ......: 8 + 2017-08-26T10:35:16: Number of spatial par's ...: 2 + 2017-08-26T10:35:16: RA .......................: 128.84 [-360,360] deg (fixed,scale=1) + 2017-08-26T10:35:16: DEC ......................: -45.18 [-90,90] deg (fixed,scale=1) + 2017-08-26T10:35:16: Number of spectral par's ..: 5 + 2017-08-26T10:35:16: Prefactor ................: 4.3553344038324e-09 +/- 4.75309654008435e-10 [1e-16,1e-06] ph/cm2/s/MeV (free,scale=1e-09,gradient) + 2017-08-26T10:35:16: Index1 ...................: -1.3435116067011 +/- 0.0311270092486095 [-0,-5] (free,scale=-1,gradient) + 2017-08-26T10:35:16: CutoffEnergy .............: 988.499237722443 +/- 181.86709257978 [1,1000000] MeV (free,scale=1000,gradient) + 2017-08-26T10:35:16: PivotEnergy ..............: 1000 [10,1000000] MeV (fixed,scale=1000,gradient) + 2017-08-26T10:35:16: Index2 ...................: 0.587129723482927 +/- 0.0279303486200403 [0.1,5] (free,scale=1,gradient) + 2017-08-26T10:35:16: Number of temporal par's ..: 1 + 2017-08-26T10:35:16: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2017-08-26T10:35:16: === GModelSky === + 2017-08-26T10:35:16: Name ......................: Galactic_diffuse + 2017-08-26T10:35:16: Instruments ...............: all + 2017-08-26T10:35:16: Instrument scale factors ..: unity + 2017-08-26T10:35:16: Observation identifiers ...: all + 2017-08-26T10:35:16: Model type ................: DiffuseSource + 2017-08-26T10:35:16: Model components ..........: "DiffuseMapCube" * "Constant" * "Constant" + 2017-08-26T10:35:16: Number of parameters ......: 3 + 2017-08-26T10:35:16: Number of spatial par's ...: 1 + 2017-08-26T10:35:16: Normalization ............: 1 [0.1,10] (fixed,scale=1,gradient) + 2017-08-26T10:35:16: Number of spectral par's ..: 1 + 2017-08-26T10:35:16: Normalization ............: 1.10915296009432 +/- 0.00503367301009883 [0.1,1000] (free,scale=1,gradient) + 2017-08-26T10:35:16: Number of temporal par's ..: 1 + 2017-08-26T10:35:16: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + 2017-08-26T10:35:16: === GModelSky === + 2017-08-26T10:35:16: Name ......................: Extragalactic_diffuse + 2017-08-26T10:35:16: Instruments ...............: all + 2017-08-26T10:35:16: Instrument scale factors ..: unity + 2017-08-26T10:35:16: Observation identifiers ...: all + 2017-08-26T10:35:16: Model type ................: DiffuseSource + 2017-08-26T10:35:16: Model components ..........: "DiffuseIsotropic" * "FileFunction" * "Constant" + 2017-08-26T10:35:16: Number of parameters ......: 3 + 2017-08-26T10:35:16: Number of spatial par's ...: 1 + 2017-08-26T10:35:16: Value ....................: 1 [0,10] (fixed,scale=1,gradient) + 2017-08-26T10:35:16: Number of spectral par's ..: 1 + 2017-08-26T10:35:16: Normalization ............: 1 [0,1000] (fixed,scale=1,gradient) + 2017-08-26T10:35:16: Number of temporal par's ..: 1 + 2017-08-26T10:35:16: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + + +The fit converged after 19 iterations with spectral parameters that are +reasonably close to those found in +`Abdo et al. 2010, ApJ, 713, 154 `_ +using the Fermi-LAT Science Tools. + + .. note:: + + The fit results can be compared to the corresponding output of gtlike + when the tool is applied to the same test data. As you will see below, + the fitted model parameters are basically identical. + + .. code-block:: bash + + Extragalactic_diffuse: + Normalization: 1 + Flux: 0.000149604 photons/cm^2/s + + Galactic_diffuse: + Value: 1.11394 +/- 0.00508393 + Flux: 0.00054415 +/- 2.48325e-06 photons/cm^2/s + + Vela: + Prefactor: 4.3344 +/- 0.65258 + Index1: 1.34634 +/- 0.0429336 + Scale: 1 + Cutoff: 1.00992 +/- 0.25512 + Index2: 0.595018 +/- 0.0398585 + TS value: 326648 + Flux: 1.0758e-05 +/- 4.33949e-08 photons/cm^2/s + + WARNING: Fit may be bad in range [100, 149.23] (MeV) + WARNING: Fit may be bad in range [182.299, 405.972] (MeV) + WARNING: Fit may be bad in range [6690.96, 8173.66] (MeV) + WARNING: Fit may be bad in range [14900.5, 18202.4] (MeV) + WARNING: Fit may be bad in range [33182.8, 40536] (MeV) + WARNING: Fit may be bad in range [49518.7, 60491.9] (MeV) + + Total number of observed counts: 202330 + Total number of model events: 202327 + + -log(Likelihood): -215606.3555 + + Elapsed CPU time: 72.33 + + You may also compute the Test Statistic value using :ref:`ctlike` by + adding the ``tscalc="1"`` attribute to the + :ref:`model definition file `: + + .. code-block:: xml + + + + + ... + + + This results in + + .. code-block:: bash + + 2017-08-26T13:05:12: === GModelSky === + 2017-08-26T13:05:12: Name ......................: Vela + 2017-08-26T13:05:12: Instruments ...............: all + 2017-08-26T13:05:12: Test Statistic ............: 326326.556902542 + + which is also close to the value determined by gtlike. diff --git a/doc/source/users/tutorials/howto/fermi/howto_fermi_prepare.rst b/doc/source/users/tutorials/howto/fermi/howto_fermi_prepare.rst new file mode 100644 index 0000000..74bd859 --- /dev/null +++ b/doc/source/users/tutorials/howto/fermi/howto_fermi_prepare.rst @@ -0,0 +1,226 @@ +.. _howto_fermi_prepare: + +Prepare Fermi-LAT data +---------------------- + + .. admonition:: What you will learn + + You will learn how to prepare Fermi-LAT data for an analysis with + ctools. + + You will find more information about the Fermi-LAT software and data + analysis at the + `Fermi Science Support Center (FSSC) `_ + web site. + + +Prerequisits +~~~~~~~~~~~~ + +The following sections assume that you have properly installed the Fermi-LAT +Science Tools. The Fermi-LAT Science Tools can be downloaded from the +`Fermi Science Support Center (FSSC) `_. + + +Get data from the Fermi-LAT data server +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Start with requesting the data you are interested in from the +`Fermi-LAT Data server `_. +As example the figure below shows a request for data around the Vela pulsar. + +.. figure:: howto_fermi_query.jpg + :width: 600px + :align: center + + *Fermi-LAT data server interface* + +After submitting the request you will get information of how to retrieve the +data: + +.. figure:: howto_fermi_data.jpg + :width: 600px + :align: center + + *Fermi-LAT data retrieval information* + +Now you download the data using the ``wget`` tool. At the same time of +downloading the data you should also download the input files for the diffuse +galactic and extragalactic model components. In the current example the +following files were downloaded: + +.. code-block:: bash + + $ wget http://fermi.gsfc.nasa.gov/FTP/fermi/data/lat/queries/L1708251631252315B1AD08_SC00.fits + $ wget http://fermi.gsfc.nasa.gov/FTP/fermi/data/lat/queries/L1708251631252315B1AD08_PH00.fits + $ wget http://fermi.gsfc.nasa.gov/FTP/fermi/data/lat/queries/L1708251631252315B1AD08_PH01.fits + $ wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/aux/gll_iem_v06.fits + $ wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/aux/iso_P8R2_SOURCE_V6_v06.txt +.. + + .. warning:: + + In the current example data within a radius of only 10 degrees around the + source of interest were selected. This was mainly for getting a small + data set. The + `Fermi Science Support Center (FSSC) `_ + recommends to use larger regions for the Fermi-LAT data analysis. + + +Data preparation +~~~~~~~~~~~~~~~~ + +First you need to combine the event lists into a single file. At the same +time you will select the region of interest for the analysis. You do this +using the ``gtselect`` tool: + +.. code-block:: bash + + $ ls *_PH* > events.txt + $ gtselect evclass=128 evtype=3 + Input FT1 file[] @events.txt + Output FT1 file[] events_fermi.fits + RA for new search center (degrees) (0:360) [INDEF] 128.84 + Dec for new search center (degrees) (-90:90) [INDEF] -45.18 + radius of new search region (degrees) (0:180) [INDEF] 10.0 + start time (MET in s) (0:) [INDEF] + end time (MET in s) (0:) [INDEF] + lower energy limit (MeV) (0:) [30] 100.0 + upper energy limit (MeV) (0:) [300000] 300000.0 + maximum zenith angle value (degrees) (0:180) [180] 90.0 + +Now you have to select from all events those which fall into periods where +data quality is good, the telescope is configured for science, and the rocking +angle is not too large. You do this using the ``gtmktime`` tool that defines +the Good Time Intervals for your analysis: + +.. code-block:: bash + + $ gtmktime + Spacecraft data file[] L1708251631252315B1AD08_SC00.fits + Filter expression[DATA_QUAL>0 && LAT_CONFIG==1 && ABS(ROCK_ANGLE)<52] + Apply ROI-based zenith angle cut[yes] no + Event data file[] events_fermi.fits + Output event file name[] events_fermi_gti.fits + +The events are now ready to be binned into a counts cube using the +``gtbin`` tool: + +.. code-block:: bash + + $ gtbin + This is gtbin version ScienceTools-10-01-01 + Type of output file (CCUBE|CMAP|LC|PHA1|PHA2|HEALPIX) [PHA2] CCUBE + Event data file name[] events_fermi_gti.fits + Output file name[] cntmap.fits + Spacecraft data file name[NONE] L1708251631252315B1AD08_SC00.fits + Size of the X axis in pixels[] 60 + Size of the Y axis in pixels[] 60 + Image scale (in degrees/pixel)[] 0.2 + Coordinate system (CEL - celestial, GAL -galactic) (CEL|GAL) [CEL] + First coordinate of image center in degrees (RA or galactic l)[] 128.84 + Second coordinate of image center in degrees (DEC or galactic b)[] -45.18 + Rotation angle of image axis, in degrees[0.] + Projection method e.g. AIT|ARC|CAR|GLS|MER|NCP|SIN|STG|TAN:[AIT] TAN + Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] + Start value for first energy bin in MeV[30] 100.0 + Stop value for last energy bin in MeV[200000] 300000.0 + Number of logarithmically uniform energy bins[] 40 + +As next step you need to compute a livetime cube which is a computational +intensive task. You do this using the ``gtltcube`` tool: + +.. code-block:: bash + + $ gtltcube zmax=90 + Event data file[] events_fermi_gti.fits + Spacecraft data file[] L1708251631252315B1AD08_SC00.fits + Output file[expCube.fits] ltcube.fits + Step size in cos(theta) (0.:1.) [0.025] + Pixel size (degrees)[1] + Working on file L1708251631252315B1AD08_SC00.fits + .....................! + +The little dots indicate the progress while the tool is computing. Once +finished you need to compute the binned exposure map using the ``gtexpcube2`` +tool: + +.. code-block:: bash + + $ gtexpcube2 + Livetime cube file[] ltcube.fits + Counts map file[] none + Output file name[] expmap.fits + Response functions to use[CALDB] P8R2_SOURCE_V6 + Size of the X axis in pixels[INDEF] 200 + Size of the Y axis in pixels[INDEF] 200 + Image scale (in degrees/pixel)[INDEF] 0.2 + Coordinate system (CEL - celestial, GAL -galactic) (CEL|GAL) [GAL] CEL + First coordinate of image center in degrees (RA or galactic l)[INDEF] 128.84 + Second coordinate of image center in degrees (DEC or galactic b)[INDEF] -45.18 + Rotation angle of image axis, in degrees[0] + Projection method e.g. AIT|ARC|CAR|GLS|MER|NCP|SIN|STG|TAN[CAR] TAN + Start energy (MeV) of first bin[INDEF] 100.0 + Stop energy (MeV) of last bin[INDEF] 300000.0 + Number of logarithmically-spaced energy bins[INDEF] 40 + Computing binned exposure map....................! + + +Generate source maps +~~~~~~~~~~~~~~~~~~~~ + +Now you are ready to produce source maps for all diffuse model components. +Start with putting the diffuse model components into a +:ref:`model definition XML file `: + +.. code-block:: xml + + + + + + + + + + + + + + + + + + + + +.. + + .. warning:: + + While ctools also implements the syntax used for Fermi-LAT Science Tools + :ref:`model definition files `, + the Fermi-LAT Science Tools naming conventions are not very homogenous, + hence we usually advocate to use the more coherent ctools syntax. For the + Fermi-LAT Science Tools you have to use however the Fermi-LAT syntax. + +Now you can generate the source maps using ``gtsrcmaps``: + +.. code-block:: bash + + $ gtsrcmaps ptsrc=no + Exposure hypercube file[] ltcube.fits + Counts map file[] cntmap.fits + Source model file[] diffuse.xml + Binned exposure map[none] expmap.fits + Source maps output file[] srcmaps.fits + Response functions[CALDB] P8R2_SOURCE_V6 + Generating SourceMap for Extragalactic_diffuse....................! + Generating SourceMap for Galactic_diffuse....................! + +And you are done with the preparation of the Fermi-LAT data. You now have +the following files in your working directory: + +* ``srcmaps.fits`` - Counts cube including source maps +* ``ltcube.fits`` - Livetime cube +* ``expmap.fits`` - Exposure map diff --git a/doc/source/users/tutorials/howto/fermi/howto_fermi_query.jpg b/doc/source/users/tutorials/howto/fermi/howto_fermi_query.jpg new file mode 100644 index 0000000..077047d --- /dev/null +++ b/doc/source/users/tutorials/howto/fermi/howto_fermi_query.jpg diff --git a/doc/source/users/tutorials/howto/fermi/howto_fermi_spectrum.png b/doc/source/users/tutorials/howto/fermi/howto_fermi_spectrum.png new file mode 100644 index 0000000..5377e79 --- /dev/null +++ b/doc/source/users/tutorials/howto/fermi/howto_fermi_spectrum.png diff --git a/doc/source/users/tutorials/howto/fermi/howto_fermi_spectrum.rst b/doc/source/users/tutorials/howto/fermi/howto_fermi_spectrum.rst new file mode 100644 index 0000000..c3ae036 --- /dev/null +++ b/doc/source/users/tutorials/howto/fermi/howto_fermi_spectrum.rst @@ -0,0 +1,83 @@ +.. _howto_fermi_spectrum: + +Generate a spectral energy distribution for a source +---------------------------------------------------- + + .. admonition:: What you will learn + + You will learn how to use :ref:`ctlike` to derive the spectral energy + density of a source from Fermi-LAT data. + +So far the :ref:`csspec` script does not work for Fermi-LAT data, but you +will learn here a method that you can use to derive the spectrum of any +source using Fermi-LAT data. For this you first need a little Python code +to transform the spectral component of the Vela source in the +:ref:`model definition file ` +from a super exponentially cut-off power law into a node function. +You do this as follows: + +.. code-block:: python + + >>> import gammalib + >>> models = gammalib.GModels('vela_results.xml') + >>> energies = gammalib.GEnergies(10, gammalib.GEnergy(100.0, 'MeV'), gammalib.GEnergy(100.0, 'GeV')) + >>> nodes = gammalib.GModelSpectralNodes() + >>> for energy in energies: + >>> nodes.append(energy, models['Vela'].spectral().eval(energy)) + >>> for par in nodes: + >>> if 'Energy' in par.name(): + >>> par.fix() + >>> if 'Intensity' in par.name(): + >>> par.min(1.0e-30) + >>> models['Vela'].spectral(nodes) + >>> models.save('model_nodes.xml') + +This will create a +:ref:`model definition file ` +with the following structure: + +.. code-block:: xml + + + + + + + + + + + + + + ... + + +Instead of the super exponentially cut-off power law the model now contains +a node function with ten intensity nodes between 100 MeV and 100 GeV, and the +:ref:`ctlike` tool will now adjust the intensities of each of the nodes which +will result in a spectrum. To do this type now: + +.. code-block:: bash + + $ ctlike debug=yes + Input event list, counts cube or observation definition XML file [obs.xml] + Input model definition XML file [models.xml] model_nodes.xml + Output model definition XML file [vela_results.xml] vela_nodes_results.xml + +This will produce the +:ref:`model definition file ` +``vela_nodes_results.xml`` with the intensities of the fitted source at +each energy. The results are displayed graphically below. + +.. figure:: howto_fermi_spectrum.png + :width: 600px + :align: center + + *Spectrum of the Vela pulsar derived using a node function* +.. + + .. note:: + + The ``NodeFunction`` spectral models implements a piece-wise power law + that is defined by the intensity values at each node. diff --git a/doc/source/users/tutorials/howto/fermi/howto_fermi_ulimit.rst b/doc/source/users/tutorials/howto/fermi/howto_fermi_ulimit.rst new file mode 100644 index 0000000..1c90ff9 --- /dev/null +++ b/doc/source/users/tutorials/howto/fermi/howto_fermi_ulimit.rst @@ -0,0 +1,61 @@ +.. _howto_fermi_ulimit: + +Derive upper flux limits for a source +------------------------------------- + + .. admonition:: What you will learn + + You will learn how to use :ref:`ctulimit` to determine upper flux limits + from Fermi-LAT data using a likelihood profile method. + +Suppose there is another potential source in the field of view but you cannot +see this source in a counts map, hence you want to derive an upper flux limit +for the source. To do this you first have to add an additional source to the +:ref:`model definition file `. An example for a source at +:math:`\alpha=122^\circ` and :math:`\delta=-42^\circ` is shown below: + +.. code-block:: xml + + + + + + + + + + + + + + + + + ... + + +Now you can run the :ref:`ctulimit` tool. Since the tool is normally tuned +for Cherenkov Telescope data analysis, energies are by default given in TeV. +To specify a reference energy of the differential flux limit of 1 GeV you +need to provide the hidden parameter ``eref=0.001``. You can also set the +flux interval for the intergal flux computation to 100 MeV - 100 GeV by +specifying ``emin=0.0001 emax=0.1`` (all values are given in TeV): + +.. code-block:: bash + + $ ctulimit eref=0.001 emin=0.0001 emax=0.1 + Input event list, counts cube or observation definition XML file [obs.xml] + Source of interest [Test] + Input model definition XML file [model_lat_ulimit.xml] + +The results of the upper limit computation can then be extracted from the +``ctulimit.log`` file that is created by :ref:`ctulimit`: + +.. code-block:: bash + + 2017-08-26T14:37:48: +=====================+ + 2017-08-26T14:37:48: | Upper limit results | + 2017-08-26T14:37:48: +=====================+ + 2017-08-26T14:37:48: Differential flux limit ...: 1.00100099727198e-13 ph/cm2/s/MeV at 0.001 TeV + 2017-08-26T14:37:48: Integral flux limit .......: 9.9999999627471e-10 ph/cm2/s within [0.0001-0.1] TeV + 2017-08-26T14:37:48: Energy flux limit .........: 1.10785216561628e-12 erg/cm2/s within [0.0001-0.1] TeV diff --git a/doc/source/users/tutorials/howto/fermi/index.rst b/doc/source/users/tutorials/howto/fermi/index.rst new file mode 100644 index 0000000..a2e1359 --- /dev/null +++ b/doc/source/users/tutorials/howto/fermi/index.rst @@ -0,0 +1,37 @@ +.. _howto_fermi: + +Analysing Fermi-LAT data +------------------------ + + .. admonition:: What you will learn + + You will learn how to analyse data from the LAT telescope about NASA's + `Fermi `_ satellite with ctools. + +You can use ctools to analyse data from the LAT telescope about NASA's +`Fermi `_ satellite. +To do this you need to be familar with the Fermi-LAT Science Tools that you +will use to prepare the Fermi-LAT data for the analysis. You can download the +Fermi-LAT Science Tools from the +`Fermi Science Support Center (FSSC) `_ +which provides also more detailed information about how to prepare the data. +What you need for a given observation are the following files: + +* ``srcmaps.fits`` - source maps file that include a counts cube and maps of the + diffuse model components +* ``ltcube.fits`` - livetime cube file that corresponds to the Good Time Intervals + of the data +* ``expmap.fits`` - exposure map file + +To demonstrate how the Fermi-LAT analysis works with ctools, here is a short +tutorial that is based on a simple analysis of the Vela pulsar. + +.. toctree:: + :maxdepth: 1 + + howto_fermi_prepare + howto_fermi_fitting + howto_fermi_butterfly + howto_fermi_spectrum + howto_fermi_ulimit + howto_fermi_errors diff --git a/doc/source/users/tutorials/iact_analysis/analysis.rst b/doc/source/users/tutorials/howto/iact/analysis.rst index 3438468..3438468 100644 --- a/doc/source/users/tutorials/iact_analysis/analysis.rst +++ b/doc/source/users/tutorials/howto/iact/analysis.rst diff --git a/doc/source/users/tutorials/iact_analysis/copy.rst b/doc/source/users/tutorials/howto/iact/copy.rst index 5ac7365..5ac7365 100644 --- a/doc/source/users/tutorials/iact_analysis/copy.rst +++ b/doc/source/users/tutorials/howto/iact/copy.rst diff --git a/doc/source/users/tutorials/iact_analysis/high_level.rst b/doc/source/users/tutorials/howto/iact/high_level.rst index c4d6d31..c4d6d31 100644 --- a/doc/source/users/tutorials/iact_analysis/high_level.rst +++ b/doc/source/users/tutorials/howto/iact/high_level.rst diff --git a/doc/source/users/tutorials/iact_analysis/index.rst b/doc/source/users/tutorials/howto/iact/index.rst index a06cccf..a06cccf 100644 --- a/doc/source/users/tutorials/iact_analysis/index.rst +++ b/doc/source/users/tutorials/howto/iact/index.rst diff --git a/doc/source/users/tutorials/iact_analysis/simulation.rst b/doc/source/users/tutorials/howto/iact/simulation.rst index 9e2ac45..9e2ac45 100644 --- a/doc/source/users/tutorials/iact_analysis/simulation.rst +++ b/doc/source/users/tutorials/howto/iact/simulation.rst diff --git a/doc/source/users/tutorials/iact_analysis/tips.rst b/doc/source/users/tutorials/howto/iact/tips.rst index 1f09b91..1f09b91 100644 --- a/doc/source/users/tutorials/iact_analysis/tips.rst +++ b/doc/source/users/tutorials/howto/iact/tips.rst diff --git a/doc/source/users/tutorials/howto/index.rst b/doc/source/users/tutorials/howto/index.rst index 0af03fb..c6d5ac7 100644 --- a/doc/source/users/tutorials/howto/index.rst +++ b/doc/source/users/tutorials/howto/index.rst @@ -21,4 +21,6 @@ These examples can be easily adapted to any data you have at your disposal. howto_lightcurve howto_phasecurve howto_exclude - howto_vo \ No newline at end of file + howto_vo + How to analyse Fermi-LAT data with ctools? + How to analyse H.E.S.S. data with ctools? \ No newline at end of file diff --git a/doc/source/users/tutorials/index.rst b/doc/source/users/tutorials/index.rst index e21ba68..015ce72 100644 --- a/doc/source/users/tutorials/index.rst +++ b/doc/source/users/tutorials/index.rst @@ -8,6 +8,5 @@ Tutorials quickstart/index beyond/index - 1dc/index howto/index - iact_analysis/index + 1dc/index \ No newline at end of file