diff --git a/notebooks/post_DL3_analysis.ipynb b/notebooks/post_DL3_analysis.ipynb index 05d73d0e76..5d663f55fe 100644 --- a/notebooks/post_DL3_analysis.ipynb +++ b/notebooks/post_DL3_analysis.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "fe60fa59", + "id": "0", "metadata": { "tags": [] }, @@ -53,7 +53,7 @@ { "cell_type": "code", "execution_count": null, - "id": "0e80c1fe", + "id": "1", "metadata": {}, "outputs": [], "source": [ @@ -70,7 +70,7 @@ { "cell_type": "code", "execution_count": null, - "id": "f1babec1", + "id": "2", "metadata": {}, "outputs": [], "source": [ @@ -109,7 +109,7 @@ { "cell_type": "code", "execution_count": null, - "id": "3a4024b6-106a-4a7b-b2f8-10f7a3981725", + "id": "3", "metadata": {}, "outputs": [], "source": [ @@ -118,7 +118,7 @@ }, { "cell_type": "markdown", - "id": "2f59f7bb-3b9c-49a0-acc7-ea3b3e0aceca", + "id": "4", "metadata": {}, "source": [ "# Input and output settings (SET BY THE USER)\n", @@ -133,7 +133,7 @@ { "cell_type": "code", "execution_count": null, - "id": "0d4ed5a5-101a-4f66-ace9-858e7c1adf83", + "id": "5", "metadata": {}, "outputs": [], "source": [ @@ -153,7 +153,7 @@ }, { "cell_type": "markdown", - "id": "c9c1308b", + "id": "6", "metadata": {}, "source": [ "# 0. Inspect the DL3 file content" @@ -162,7 +162,7 @@ { "cell_type": "code", "execution_count": null, - "id": "ae7c573e-71da-4984-a3c9-cbda09a832f8", + "id": "7", "metadata": {}, "outputs": [], "source": [ @@ -174,7 +174,7 @@ { "cell_type": "code", "execution_count": null, - "id": "242c39c3-764e-45de-b7a8-db63e5b815bc", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -185,7 +185,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e3ac3e87-5d07-40e1-9b71-3c5e04ada5e4", + "id": "9", "metadata": {}, "outputs": [], "source": [ @@ -196,7 +196,7 @@ { "cell_type": "code", "execution_count": null, - "id": "ff6d1ab8-a995-44d6-8cd2-42a58887b946", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -208,7 +208,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e45b1cee-bfb6-4e7e-8f9d-cf1effb7f245", + "id": "11", "metadata": {}, "outputs": [], "source": [ @@ -220,7 +220,7 @@ }, { "cell_type": "markdown", - "id": "cbc4490a-23bd-45ad-81cc-8a8f9430bedb", + "id": "12", "metadata": {}, "source": [ "# 1. Read the DL3 index files and load the data" @@ -228,7 +228,7 @@ }, { "cell_type": "markdown", - "id": "658cbcf5", + "id": "13", "metadata": {}, "source": [ "Let's load now all the DL3 files together.\n", @@ -240,7 +240,7 @@ { "cell_type": "code", "execution_count": null, - "id": "08dc4bdf", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -250,7 +250,7 @@ { "cell_type": "code", "execution_count": null, - "id": "29b54009-d0cd-4c06-a68b-98ebdaa9b6ae", + "id": "15", "metadata": { "tags": [] }, @@ -264,9 +264,8 @@ { "cell_type": "code", "execution_count": null, - "id": "aa584be6", + "id": "16", "metadata": { - "scrolled": true, "tags": [] }, "outputs": [], @@ -278,7 +277,7 @@ }, { "cell_type": "markdown", - "id": "f683b5ad", + "id": "17", "metadata": {}, "source": [ "# 2. Selection filters for the observations\n", @@ -291,7 +290,7 @@ { "cell_type": "code", "execution_count": null, - "id": "40149928", + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -308,7 +307,7 @@ { "cell_type": "code", "execution_count": null, - "id": "a2527c35", + "id": "19", "metadata": {}, "outputs": [], "source": [ @@ -323,9 +322,8 @@ { "cell_type": "code", "execution_count": null, - "id": "49f8226c", + "id": "20", "metadata": { - "scrolled": true, "tags": [] }, "outputs": [], @@ -346,7 +344,7 @@ { "cell_type": "code", "execution_count": null, - "id": "37d39c99", + "id": "21", "metadata": {}, "outputs": [], "source": [ @@ -357,7 +355,7 @@ }, { "cell_type": "markdown", - "id": "5a4220a9", + "id": "22", "metadata": {}, "source": [ "# 3. Define Target position and energy ranges for reconstructed events" @@ -366,7 +364,7 @@ { "cell_type": "code", "execution_count": null, - "id": "86fa0999", + "id": "23", "metadata": {}, "outputs": [], "source": [ @@ -376,7 +374,7 @@ }, { "cell_type": "markdown", - "id": "6d00c086-8be2-4815-a450-a4054d5fa41f", + "id": "24", "metadata": {}, "source": [ "If we were using global theta cut, this is the way of getting the theta cut (`RAD_MAX` key) used for the IRF production.\n", @@ -386,7 +384,7 @@ }, { "cell_type": "markdown", - "id": "7cc29858-9b8e-49f1-9e7a-daffe127091f", + "id": "25", "metadata": {}, "source": [ "```\n", @@ -403,7 +401,7 @@ { "cell_type": "code", "execution_count": null, - "id": "bc0afac4-476b-4708-b9ef-8b19f807c371", + "id": "26", "metadata": {}, "outputs": [], "source": [ @@ -414,7 +412,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e3f9d096", + "id": "27", "metadata": {}, "outputs": [], "source": [ @@ -468,7 +466,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e34f4149", + "id": "28", "metadata": {}, "outputs": [], "source": [ @@ -478,7 +476,7 @@ }, { "cell_type": "markdown", - "id": "1ecb339f", + "id": "29", "metadata": {}, "source": [ "# 4. Define the base Map geometries for creating the SpectrumDataset\n", @@ -489,7 +487,7 @@ { "cell_type": "code", "execution_count": null, - "id": "2ac5aee7", + "id": "30", "metadata": {}, "outputs": [], "source": [ @@ -504,7 +502,7 @@ }, { "cell_type": "markdown", - "id": "06651bb8", + "id": "31", "metadata": {}, "source": [ "# 5. Data Reduction chain\n", @@ -527,7 +525,7 @@ { "cell_type": "code", "execution_count": null, - "id": "92b8885e", + "id": "32", "metadata": {}, "outputs": [], "source": [ @@ -546,7 +544,7 @@ }, { "cell_type": "markdown", - "id": "0a6bb7be-17a2-4e75-82dd-e2ebf89ce62b", + "id": "33", "metadata": {}, "source": [ "The following makers can be tuned to check the final Dataset to be used." @@ -555,7 +553,7 @@ { "cell_type": "code", "execution_count": null, - "id": "83510236-5380-422d-96c5-dc09a3134d31", + "id": "34", "metadata": {}, "outputs": [], "source": [ @@ -566,7 +564,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9c49a8f0-ca04-4df4-8a5e-a34d32fd51c6", + "id": "35", "metadata": {}, "outputs": [], "source": [ @@ -581,7 +579,7 @@ { "cell_type": "code", "execution_count": null, - "id": "ee64b572-14b9-4729-8b59-9a3a2cc17381", + "id": "36", "metadata": {}, "outputs": [], "source": [ @@ -592,7 +590,7 @@ }, { "cell_type": "markdown", - "id": "30fc1f88-5833-4ff4-bbe4-1a42b2312aa8", + "id": "37", "metadata": {}, "source": [ "For other arguments and options, check the SafeMaskMaker documentation:\n", @@ -601,7 +599,7 @@ }, { "cell_type": "markdown", - "id": "4c37d609", + "id": "38", "metadata": {}, "source": [ "# 6. Generate the Spectrum Dataset for all observations" @@ -610,10 +608,8 @@ { "cell_type": "code", "execution_count": null, - "id": "48b17e0c", - "metadata": { - "scrolled": true - }, + "id": "39", + "metadata": {}, "outputs": [], "source": [ "%%time\n", @@ -648,7 +644,7 @@ }, { "cell_type": "markdown", - "id": "899f6974-cb53-492b-9655-c0083a2ee7d4", + "id": "40", "metadata": {}, "source": [ "Have a look at the information from the first dataset" @@ -657,7 +653,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d9b75eea", + "id": "41", "metadata": {}, "outputs": [], "source": [ @@ -666,7 +662,7 @@ }, { "cell_type": "markdown", - "id": "7158fb21-8231-480f-92e0-e6163ed69c6e", + "id": "42", "metadata": {}, "source": [ "## 6.1 Stack datasets" @@ -674,7 +670,7 @@ }, { "cell_type": "markdown", - "id": "780be7a1-b799-4d51-af07-99f075e93978", + "id": "43", "metadata": {}, "source": [ "At this point, all the datasets can be stacked in a single dataset (SpectrumDatasetOnOff object), which is usually faster to handle.\n", @@ -685,7 +681,7 @@ { "cell_type": "code", "execution_count": null, - "id": "0bea45f4-de4e-4af3-b141-d3ebb01953f4", + "id": "44", "metadata": {}, "outputs": [], "source": [ @@ -694,7 +690,7 @@ }, { "cell_type": "markdown", - "id": "710624fb-c26d-495d-8712-16d9a57d8264", + "id": "45", "metadata": {}, "source": [ "Alternatively, you can continue directly using the individual `Datasets` as produced in the above loop. Then you would perform a joint likelihood fit to all observations." @@ -702,7 +698,7 @@ }, { "cell_type": "markdown", - "id": "fea40a2a", + "id": "46", "metadata": {}, "source": [ "# 7. Some plots with the given Dataset" @@ -710,7 +706,7 @@ }, { "cell_type": "markdown", - "id": "c863d7b9-b72e-4717-8a70-0da09434d0b5", + "id": "47", "metadata": {}, "source": [ "Plot the target position and center of the OFF regions used for the calculation of the gamma-ray excess." @@ -719,7 +715,7 @@ { "cell_type": "code", "execution_count": null, - "id": "dcaea147", + "id": "48", "metadata": {}, "outputs": [], "source": [ @@ -731,7 +727,7 @@ }, { "cell_type": "raw", - "id": "3bf4ab3f", + "id": "49", "metadata": {}, "source": [ "# For source-dependent analysis, check the reconstructed position of all the events, \n", @@ -748,10 +744,8 @@ { "cell_type": "code", "execution_count": null, - "id": "5bac0518", - "metadata": { - "scrolled": true - }, + "id": "50", + "metadata": {}, "outputs": [], "source": [ "info_table = datasets.info_table(cumulative=True)\n", @@ -761,7 +755,7 @@ { "cell_type": "code", "execution_count": null, - "id": "899bc210", + "id": "51", "metadata": {}, "outputs": [], "source": [ @@ -794,10 +788,8 @@ { "cell_type": "code", "execution_count": null, - "id": "4f9fdb9c", - "metadata": { - "scrolled": true - }, + "id": "52", + "metadata": {}, "outputs": [], "source": [ "%%time\n", @@ -828,7 +820,7 @@ }, { "cell_type": "markdown", - "id": "d8e5f672", + "id": "53", "metadata": {}, "source": [ "# 8. Write all datasets into OGIP files\n", @@ -839,10 +831,8 @@ { "cell_type": "code", "execution_count": null, - "id": "1177b138", - "metadata": { - "scrolled": true - }, + "id": "54", + "metadata": {}, "outputs": [], "source": [ "for dataset in datasets:\n", @@ -854,7 +844,7 @@ { "cell_type": "code", "execution_count": null, - "id": "aa00c214", + "id": "55", "metadata": {}, "outputs": [], "source": [ @@ -868,7 +858,7 @@ }, { "cell_type": "markdown", - "id": "57330fc9", + "id": "56", "metadata": {}, "source": [ "# 9. Get Pivot energy to fix the reference energy and define the Spectrum Model" @@ -876,7 +866,7 @@ }, { "cell_type": "markdown", - "id": "d46242bf-1a66-4dfe-9e42-9b6aa51033ca", + "id": "57", "metadata": {}, "source": [ "As of Gammapy v1.2 it is possible to calculate the pivot energy for any spectral model not just for the power law like it is show below" @@ -885,7 +875,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8f8e421f", + "id": "58", "metadata": {}, "outputs": [], "source": [ @@ -920,7 +910,7 @@ { "cell_type": "code", "execution_count": null, - "id": "84a8072b", + "id": "59", "metadata": {}, "outputs": [], "source": [ @@ -932,7 +922,7 @@ { "cell_type": "code", "execution_count": null, - "id": "ea893639", + "id": "60", "metadata": {}, "outputs": [], "source": [ @@ -963,7 +953,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e31a06ac", + "id": "61", "metadata": {}, "outputs": [], "source": [ @@ -973,7 +963,7 @@ }, { "cell_type": "markdown", - "id": "69f4de45", + "id": "62", "metadata": {}, "source": [ "# 10. Spectral Fitting\n", @@ -985,7 +975,7 @@ { "cell_type": "code", "execution_count": null, - "id": "3f74af6f", + "id": "63", "metadata": {}, "outputs": [], "source": [ @@ -996,7 +986,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4ca480c1", + "id": "64", "metadata": {}, "outputs": [], "source": [ @@ -1006,10 +996,18 @@ "model_best = model_lp" ] }, + { + "cell_type": "markdown", + "id": "65", + "metadata": {}, + "source": [ + "Have a look at the fitting results:" + ] + }, { "cell_type": "code", "execution_count": null, - "id": "c3aa71d1", + "id": "66", "metadata": {}, "outputs": [], "source": [ @@ -1019,49 +1017,100 @@ { "cell_type": "code", "execution_count": null, - "id": "16bc612a", + "id": "67", "metadata": {}, "outputs": [], "source": [ - "# Compute the Flux Points after Fitting the model\n", - "# We do not do too many optimizations here. \n", - "# If one wants, can try and check the various attributes of the Estimator\n", - "fpe = FluxPointsEstimator(\n", - " energy_edges=energy_fit_edges, \n", - " reoptimize = False, # Re-optimizing other free model parameters (not belonging to the source)\n", - " source=obj_name,\n", - " selection_optional=\"all\" # Estimates asymmetric errors, upper limits and fit statistic profiles\n", - ")\n", - "flux_points = fpe.run(datasets=stacked_dataset)\n", - "\n", - "flux_points_dataset = FluxPointsDataset(\n", - " data=flux_points, models=model_best\n", - ")" + "model_best.to_dict()['spectral']['parameters']" ] }, { "cell_type": "code", "execution_count": null, - "id": "7b85bdb3", + "id": "68", "metadata": {}, "outputs": [], "source": [ "result" ] }, + { + "cell_type": "markdown", + "id": "69", + "metadata": {}, + "source": [ + "## Compute the Flux Points after Fitting the model\n", + "\n", + "See the various attributes of the `Estimator` in Gammapy documentation (https://docs.gammapy.org/1.3/tutorials/api/estimators.html)." + ] + }, { "cell_type": "code", "execution_count": null, - "id": "08292402", + "id": "70", "metadata": {}, "outputs": [], "source": [ - "model_best.to_dict()['spectral']['parameters']" + "fpe = FluxPointsEstimator(\n", + " energy_edges=energy_fit_edges, \n", + " reoptimize = False, # Re-optimizing other free model parameters (not belonging to the source)\n", + " source=obj_name,\n", + " selection_optional=\"all\" # Estimates asymmetric errors, upper limits and fit statistic profiles\n", + ")\n", + "\n", + "flux_points = fpe.run(datasets=stacked_dataset)\n", + "\n", + "flux_points_dataset = FluxPointsDataset(\n", + " data=flux_points, models=model_best\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "71", + "metadata": {}, + "source": [ + "### Calculating upper limits for dim sources\n", + "\n", + "To avoid NaN values when calculating flux point upper limits for dim sources, you may need to widen the default\n", + "normalization range that `FluxPointsEstimator` and `LightCurveEstimator` can scan (by default set between 0.2 and 5).\n", + "\n", + "In Gammapy 1.3 you can change it like this:\n", + "\n", + "```python\n", + " norm = Parameter(name=\"norm\", value=1.0)\n", + "```\n", + "\n", + "Include in the `FluxPointsEstimator` definition as one of its arguments and then set the value of the parameter:\n", + "\n", + "```python\n", + "fpe = FluxPointsEstimator(\n", + " ...,\n", + " norm=norm,\n", + ")\n", + "fpe.norm.scan_min = 0.05\n", + "fpe.norm.scan_max = 20\n", + "```\n", + "\n", + "Finally run the estimator:\n", + "\n", + "```python\n", + "flux_points = fpe.run(datasets=datasets)\n", + "```\n", + "\n", + "In older versions of Gammapy (