Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

New plots for calibration benchmarking #43

Merged
merged 3 commits into from
Apr 6, 2020
Merged

New plots for calibration benchmarking #43

Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
597f29a
Added plots for pedestals and dc_to_phe conversion factor.
Feb 11, 2020
7758578
Added information about reference simtel file and missing images.
Feb 17, 2020
0b0412f
Simplified function 'load_reset_simtel'.
Feb 25, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
231 changes: 223 additions & 8 deletions notebooks/DL1/calibration.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,11 +32,11 @@
"\n",
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the purpose of this function? Wouldn't it just work to pass max_events to it, and if it's none, it will work fine (no need for the if-statement, I think)


Reply via ReviewNB

"For the moment this notebook is optimized to work only on files produced from LSTCam + NectarCam telescope configurations. \n",
"As with any other part of _protopipe_ and being part of the official repository, this notebook can be further developed by any interested contributor. \n",
"The execution of this notebook is not currently automatic, it must be done locally by the user - preferably _before_ pushing a pull-request.\n",
"The execution of this notebook is not currently automatic, it must be done locally by the user - preferably _before_ pushing a pull-request. \n",
"**IMPORTANT:** Please, if you wish to contribute to this notebook, before pushing anything to your branch (better even before opening the PR) clear all the output and remove any local directory paths that you used for testing (leave empty strings). The file used shouud always be _gamma_20deg_180deg_run100___cta-prod3-demo-2147m-LaPalma-baseline.simtel.gz_ until Prod5.\n",
"\n",
"**TODO:** \n",
"* R1-data level --> pedestal vs pix_id\n",
"* R1-data level --> dc_to_phe vs pix_id"
"* ..."
]
},
{
Expand All @@ -63,7 +63,10 @@
"import matplotlib.pyplot as plt\n",
"import matplotlib.colors as colors\n",
"from matplotlib.colors import LogNorm\n",
"import matplotlib.ticker as ticker"
"import matplotlib.ticker as ticker\n",
"\n",
"from ctapipe.io import event_source\n",
"from ctapipe.instrument import CameraGeometry"
]
},
{
Expand All @@ -73,6 +76,43 @@
"## Functions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Add statistical information to a plot"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def add_stats(data, ax, x = 0.70, y = 0.85, color = \"black\"):\n",
" \"\"\"Add a textbox containing statistical information.\"\"\"\n",
" mu = data.mean()\n",
" median = np.median(data)\n",
" sigma = data.std()\n",
" textstr = '\\n'.join((\n",
" r'$\\mu=%.2f$' % (mu, ),\n",
" r'$\\mathrm{median}=%.2f$' % (median, ),\n",
" r'$\\sigma=%.2f$' % (sigma, )))\n",
"\n",
" # these are matplotlib.patch.Patch properties\n",
" props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)\n",
"\n",
" # place a text box in upper left in axes coords\n",
" ax.text(x, y, \n",
" textstr, \n",
" transform=ax.transAxes, \n",
" fontsize=10,\n",
" horizontalalignment='left',\n",
" verticalalignment='center', \n",
" bbox=props,\n",
" color=color)"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down Expand Up @@ -111,6 +151,15 @@
"metadata": {},
"outputs": [],
"source": [
"def load_reset_simtel(indir = \"./\", fileName = \"gamma_20deg_180deg_run100___cta-prod3-demo-2147m-LaPalma-baseline.simtel.gz\", max_events=None, config=\"test\"):\n",
" \"\"\"(Re)load the simtel file for all events and telescopes.\"\"\"\n",
" if max_events == None:\n",
" source = event_source(input_url=f\"{indir}/{fileName}\")\n",
" else:\n",
" source = event_source(input_url=f\"{indir}/{fileName}\", max_events=max_events)\n",
" suffix = config # all generated plots will have this as a suffix in their name\n",
" return source, suffix\n",
"\n",
"def load_reset_images(indir = \"./\", fileName = \"images.h5\", config=\"test\"):\n",
" \"\"\"(Re)load the file containing the images and extract the data per telescope type.\"\"\"\n",
" # load DL1 images\n",
Expand Down Expand Up @@ -152,21 +201,80 @@
"### Load"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### MonteCarlo data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# fill with the correct path, filename of the generated file in your system\n",
"data_LST, data_MST, config = load_reset_images()"
"# load every time you want to plot simtel-related information....\n",
"indir = \"\"\n",
"infile = \"gamma_20deg_180deg_run100___cta-prod3-demo-2147m-LaPalma-baseline.simtel.gz\"\n",
"source, config = load_reset_simtel(indir=indir,\n",
" fileName=infile,\n",
" max_events=2, # 2nd event is 1st to trigger both cameras\n",
" config=\"test\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Setup"
"#### Pipeline-processed data up to DL1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**IMPORTANT Basic information about the reference simtel file** \n",
"The file used in these benchmarks is \n",
"_gamma_20deg_180deg_run100___cta-prod3-demo-2147m-LaPalma-baseline.simtel.gz_ \n",
"and has the following basic features when NO selection is applied,\n",
"* number of simulated showers = 9793\n",
"* number of images (LST + MST) = 44401\n",
"* min number of triggered telescopes per shower = 2"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# fill with the your path, filename of the generated file in your system + config name\n",
"data_LST, data_MST, config = load_reset_images(indir=\"\",\n",
" fileName=\"\",\n",
" config=\"\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"missing_images = 44401 - len(data_LST.col(\"mc_phe_image\")) - len(data_MST.col(\"mc_phe_image\"))\n",
"if missing_images:\n",
" print(f\"WARNING: it appears you are missing {missing_images} images!\")\n",
" print(f\"This corresponds to about {missing_images*100/44401:.0f}% of the total statistics.\")\n",
" print(\"Please, check that:\")\n",
" print(\"* either you have enabled some cuts in analysis.yaml,\")\n",
" print(\"* or you are not considering some events in your analysis when you write to file.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Setup of data read from pipeline files"
]
},
{
Expand Down Expand Up @@ -228,6 +336,113 @@
"## Plots"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### R1-level information"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Pedestals"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for event in source:\n",
" triggered_telescopes = np.asarray(list(event.r0.tels_with_data))\n",
" if (triggered_telescopes > 5).any() and (triggered_telescopes < 5).any():\n",
" lst_found = 0\n",
" for tel_id in triggered_telescopes:\n",
" cam_id = event.inst.subarray.tel[tel_id].camera.cam_id\n",
" pix_ids = event.inst.subarray.tel[tel_id].camera.pix_id\n",
" pedestals = event.mc.tel[tel_id].pedestal\n",
" if (lst_found == 1) and (cam_id == \"LSTCam\"):\n",
" continue\n",
" elif (lst_found == 0) and (cam_id == \"LSTCam\"):\n",
" fig, (ax1, ax2) = plt.subplots(1,2, figsize=(12, 5), tight_layout=False, sharey=True)\n",
" ax1.set_xlabel(\"Pixel ID\")\n",
" ax1.set_ylabel(\"Pedestal ADC counts\")\n",
" p1 = ax1.plot(pix_ids, pedestals[1], label=\"High gain\")\n",
" p2 = ax1.plot(pix_ids, pedestals[0], label=\"Low gain\")\n",
" ax2.hist(pedestals[1], bins = 100, orientation=\"horizontal\", label=\"pedestals HG\")\n",
" add_stats(pedestals[1], ax2, x = 0.55, y = 0.10, color = p1[0].get_color())\n",
" ax2.hist(pedestals[0], bins = 100, orientation=\"horizontal\", label=\"pedestals LG\")\n",
" add_stats(pedestals[0], ax2, x = 0.55, y = 0.25, color = p2[0].get_color())\n",
" ax1.legend()\n",
" ax2.legend()\n",
" lst_found = 1\n",
" fig.savefig(f\"./plots_calibration/pedestalsVSpixelids_{cam_id}_protopipe_{config}.png\")\n",
" elif (lst_found == 1) and (cam_id == \"NectarCam\"):\n",
" fig, (ax1, ax2) = plt.subplots(1,2, figsize=(12, 5), tight_layout=False, sharey=True)\n",
" ax1.set_xlabel(\"Pixel ID\")\n",
" ax1.set_ylabel(\"Pedestal ADC counts\")\n",
" p1 = ax1.plot(pix_ids, pedestals[1], label=\"High gain\")\n",
" p2 = ax1.plot(pix_ids, pedestals[0], label=\"Low gain\")\n",
" ax2.hist(pedestals[1], bins = 100, orientation=\"horizontal\", label=\"pedestals HG\")\n",
" add_stats(pedestals[1], ax2, x = 0.55, y = 0.10, color = p1[0].get_color())\n",
" ax2.hist(pedestals[0], bins = 100, orientation=\"horizontal\", label=\"pedestals LG\")\n",
" add_stats(pedestals[0], ax2, x = 0.55, y = 0.25, color = p2[0].get_color())\n",
" ax1.legend()\n",
" ax2.legend()\n",
" fig.savefig(f\"./plots_calibration/pedestalsVSpixelids_{cam_id}_protopipe_{config}.png\")\n",
" break"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### DC ---> PHE"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"# no cycle over events since we already have a good event from the previous cell\n",
"triggered_telescopes = np.asarray(list(event.r0.tels_with_data))\n",
"if (triggered_telescopes > 5).any() and (triggered_telescopes < 5).any():\n",
" lst_found = 0\n",
" for tel_id in triggered_telescopes:\n",
" cam_id = event.inst.subarray.tel[tel_id].camera.cam_id\n",
" pix_ids = event.inst.subarray.tel[tel_id].camera.pix_id\n",
" dc_to_pe = event.mc.tel[tel_id].dc_to_pe\n",
" if (lst_found == 1) and (cam_id == \"LSTCam\"):\n",
" continue\n",
" elif (lst_found == 0) and (cam_id == \"LSTCam\"):\n",
" lst_found = 1\n",
" dc_to_pe_channels_lst = dc_to_pe\n",
" elif (lst_found == 1) and (cam_id == \"NectarCam\"):\n",
" dc_to_pe_channels_mst = dc_to_pe\n",
" break\n",
"# plot channel-wise\n",
"for i, gain in enumerate([\"High\", \"Low\"]):\n",
" fig, (ax1, ax2) = plt.subplots(1,2, figsize=(12, 5), tight_layout=False, sharey=True)\n",
" ax1.set_xlabel(\"Pixel ID\")\n",
" ax2.set_ylabel(f\"DC to PHE - {gain} gain\")\n",
" p1 = ax1.plot(pix_ids, dc_to_pe_channels_lst[i], label=\"LSTCam\")\n",
" p2 = ax1.plot(pix_ids, dc_to_pe_channels_mst[i], label=\"NectarCam\")\n",
" ax2.hist(dc_to_pe_channels_lst[i], bins = 100, orientation=\"horizontal\", label=\"pedestals HG\")\n",
" add_stats(dc_to_pe_channels_lst[i], ax2, x = 0.55, y = 0.30, color = p1[0].get_color())\n",
" ax2.hist(dc_to_pe_channels_mst[i], bins = 100, orientation=\"horizontal\", label=\"pedestals LG\")\n",
" add_stats(dc_to_pe_channels_mst[i], ax2, x = 0.55, y = 0.55, color = p2[0].get_color())\n",
" ax1.legend()\n",
" ax2.legend()\n",
" fig.savefig(f\"./plots_calibration/dcTophe{gain}GainVSpixelids_LSTCam+NectarCam_protopipe_{config}.png\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down Expand Up @@ -558,7 +773,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
"version": "3.7.4"
}
},
"nbformat": 4,
Expand Down