.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/tutorials/raw_data_exploration.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_tutorials_raw_data_exploration.py: Exploring Raw Data ================== Here are just some very simple examples of going through and inspecting the raw data, and making some plots using ``ctapipe``. The data explored here are *raw Monte Carlo* data, which is Data Level “R0” in CTAO terminology (e.g. it is before any processing that would happen inside a Camera or off-line) .. GENERATED FROM PYTHON SOURCE LINES 15-17 Setup: .. GENERATED FROM PYTHON SOURCE LINES 17-29 .. code-block:: Python :lineno-start: 18 import numpy as np from matplotlib import pyplot as plt from scipy import signal from ctapipe.io import EventSource from ctapipe.utils import get_dataset_path from ctapipe.visualization import CameraDisplay # %matplotlib inline .. GENERATED FROM PYTHON SOURCE LINES 30-42 To read SimTelArray format data, ctapipe uses the ``pyeventio`` library (which is installed automatically along with ctapipe). The following lines however will load any data known to ctapipe (multiple ``EventSources`` are implemented, and chosen automatically based on the type of the input file. All data access first starts with an ``EventSource``, and here we use a helper function ``event_source`` that constructs one. The resulting ``source`` object can be iterated over like a list of events. We also here use an ``EventSeeker`` which provides random-access to the source (by seeking to the given event ID or number) .. GENERATED FROM PYTHON SOURCE LINES 42-46 .. code-block:: Python :lineno-start: 43 source = EventSource(get_dataset_path("gamma_prod5.simtel.zst"), max_events=5) .. GENERATED FROM PYTHON SOURCE LINES 47-53 Explore the contents of an event -------------------------------- note that the R0 level is the raw data that comes out of a camera, and also the lowest level of monte-carlo data. .. GENERATED FROM PYTHON SOURCE LINES 53-60 .. code-block:: Python :lineno-start: 54 # so we can advance through events one-by-one event_iterator = iter(source) event = next(event_iterator) .. GENERATED FROM PYTHON SOURCE LINES 61-64 the event is just a class with a bunch of data items in it. You can see a more compact representation via: .. GENERATED FROM PYTHON SOURCE LINES 64-68 .. code-block:: Python :lineno-start: 65 event.r0 .. rst-class:: sphx-glr-script-out .. code-block:: none ctapipe.containers.R0Container: tel[*]: map of tel_id to R0CameraContainer with default None .. GENERATED FROM PYTHON SOURCE LINES 69-73 printing the event structure, will currently print the value all items under it (so you get a lot of output if you print a high-level container): .. GENERATED FROM PYTHON SOURCE LINES 73-76 .. code-block:: Python :lineno-start: 74 print(event.simulation.shower) .. rst-class:: sphx-glr-script-out .. code-block:: none {'alt': , 'az': , 'core_x': , 'core_y': , 'energy': , 'h_first_int': , 'shower_primary_id': 0, 'starting_grammage': , 'x_max': } .. GENERATED FROM PYTHON SOURCE LINES 77-80 .. code-block:: Python :lineno-start: 77 print(event.r0.tel.keys()) .. rst-class:: sphx-glr-script-out .. code-block:: none dict_keys([9, 14, 104]) .. GENERATED FROM PYTHON SOURCE LINES 81-83 note that the event has 3 telescopes in it: Let’s try the next one: .. GENERATED FROM PYTHON SOURCE LINES 83-88 .. code-block:: Python :lineno-start: 84 event = next(event_iterator) print(event.r0.tel.keys()) .. rst-class:: sphx-glr-script-out .. code-block:: none dict_keys([26, 27, 29, 47, 53, 59, 69, 73, 121, 122, 124, 148, 162, 166]) .. GENERATED FROM PYTHON SOURCE LINES 89-92 now, we have a larger event with many telescopes… Let’s look at one of them: .. GENERATED FROM PYTHON SOURCE LINES 92-98 .. code-block:: Python :lineno-start: 93 teldata = event.r0.tel[26] print(teldata) teldata .. rst-class:: sphx-glr-script-out .. code-block:: none {'waveform': array([[[267, 276, 282, ..., 237, 246, 259], [254, 249, 254, ..., 210, 208, 208], [214, 233, 254, ..., 271, 272, 272], ..., [242, 236, 228, ..., 242, 252, 259], [245, 236, 230, ..., 248, 246, 248], [229, 231, 233, ..., 236, 238, 234]]], shape=(1, 1764, 25), dtype=uint16)} ctapipe.containers.R0CameraContainer: waveform: numpy array containing ADC samples: (n_channels, n_pixels, n_samples) with default None .. GENERATED FROM PYTHON SOURCE LINES 99-103 Note that some values are unit quantities (``astropy.units.Quantity``) or angular quantities (``astropy.coordinates.Angle``), and you can easily manipulate them: .. GENERATED FROM PYTHON SOURCE LINES 103-106 .. code-block:: Python :lineno-start: 104 event.simulation.shower.energy .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 107-109 .. code-block:: Python :lineno-start: 107 event.simulation.shower.energy.to("GeV") .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 110-112 .. code-block:: Python :lineno-start: 110 event.simulation.shower.energy.to("J") .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 113-115 .. code-block:: Python :lineno-start: 113 event.simulation.shower.alt .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 116-119 .. code-block:: Python :lineno-start: 116 print("Altitude in degrees:", event.simulation.shower.alt.deg) .. rst-class:: sphx-glr-script-out .. code-block:: none Altitude in degrees: 69.99999967119774 .. GENERATED FROM PYTHON SOURCE LINES 120-129 Look for signal pixels in a camera ---------------------------------- again, ``event.r0.tel[x]`` contains a data structure for the telescope data, with some fields like ``waveform``. Let’s make a 2D plot of the sample data (sample vs pixel), so we can see if we see which pixels contain Cherenkov light signals: .. GENERATED FROM PYTHON SOURCE LINES 129-136 .. code-block:: Python :lineno-start: 130 plt.pcolormesh(teldata.waveform[0]) # note the [0] is for channel 0 plt.colorbar() plt.xlabel("sample number") plt.ylabel("Pixel_id") .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_001.png :alt: raw data exploration :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(29.472222222222214, 0.5, 'Pixel_id') .. GENERATED FROM PYTHON SOURCE LINES 137-140 Let’s zoom in to see if we can identify the pixels that have the Cherenkov signal in them .. GENERATED FROM PYTHON SOURCE LINES 140-149 .. code-block:: Python :lineno-start: 141 plt.pcolormesh(teldata.waveform[0]) plt.colorbar() plt.ylim(700, 750) plt.xlabel("sample number") plt.ylabel("pixel_id") print("waveform[0] is an array of shape (N_pix,N_slice) =", teldata.waveform[0].shape) .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_002.png :alt: raw data exploration :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none waveform[0] is an array of shape (N_pix,N_slice) = (1764, 25) .. GENERATED FROM PYTHON SOURCE LINES 150-154 Now we can really see that some pixels have a signal in them! Lets look at a 1D plot of pixel 270 in channel 0 and see the signal: .. GENERATED FROM PYTHON SOURCE LINES 154-159 .. code-block:: Python :lineno-start: 155 trace = teldata.waveform[0][719] plt.plot(trace, drawstyle="steps") .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_003.png :alt: raw data exploration :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 160-164 Great! It looks like a *standard Cherenkov signal*! Let’s take a look at several traces to see if the peaks area aligned: .. GENERATED FROM PYTHON SOURCE LINES 164-172 .. code-block:: Python :lineno-start: 165 for pix_id in range(718, 723): plt.plot( teldata.waveform[0][pix_id], label="pix {}".format(pix_id), drawstyle="steps" ) plt.legend() .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_004.png :alt: raw data exploration :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 173-184 Look at the time trace from a Camera Pixel ------------------------------------------ ``ctapipe.calib.camera`` includes classes for doing automatic trace integration with many methods, but before using that, let’s just try to do something simple! Let’s define the integration windows first: By eye, they seem to be reaonsable from sample 8 to 13 for signal, and 20 to 29 for pedestal (which we define as the sum of all noise: NSB + electronic) .. GENERATED FROM PYTHON SOURCE LINES 184-192 .. code-block:: Python :lineno-start: 185 for pix_id in range(718, 723): plt.plot(teldata.waveform[0][pix_id], "+-") plt.fill_betweenx([0, 1600], 19, 24, color="red", alpha=0.3, label="Ped window") plt.fill_betweenx([0, 1600], 5, 9, color="green", alpha=0.3, label="Signal window") plt.legend() .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_005.png :alt: raw data exploration :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 193-201 Do a very simplisitic trace analysis ------------------------------------ Now, let’s for example calculate a signal and background in a the fixed windows we defined for this single event. Note we are ignoring the fact that cameras have 2 gains, and just using a single gain (channel 0, which is the high-gain channel): .. GENERATED FROM PYTHON SOURCE LINES 201-206 .. code-block:: Python :lineno-start: 202 data = teldata.waveform[0] peds = data[:, 19:24].mean(axis=1) # mean of samples 20 to 29 for all pixels sums = data[:, 5:9].sum(axis=1) / (13 - 8) # simple sum integration .. GENERATED FROM PYTHON SOURCE LINES 207-211 .. code-block:: Python :lineno-start: 207 phist = plt.hist(peds, bins=50, range=[0, 150]) plt.title("Pedestal Distribution of all pixels for a single event") .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_006.png :alt: Pedestal Distribution of all pixels for a single event :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'Pedestal Distribution of all pixels for a single event') .. GENERATED FROM PYTHON SOURCE LINES 212-215 let’s now take a look at the pedestal-subtracted sums and a pedestal-subtracted signal: .. GENERATED FROM PYTHON SOURCE LINES 215-221 .. code-block:: Python :lineno-start: 216 plt.plot(sums - peds) plt.xlabel("pixel id") plt.ylabel("Pedestal-subtracted Signal") .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_007.png :alt: raw data exploration :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(26.722222222222214, 0.5, 'Pedestal-subtracted Signal') .. GENERATED FROM PYTHON SOURCE LINES 222-226 Now, we can clearly see that the signal is centered at 0 where there is no Cherenkov light, and we can also clearly see the shower around pixel 250. .. GENERATED FROM PYTHON SOURCE LINES 226-233 .. code-block:: Python :lineno-start: 227 # we can also subtract the pedestals from the traces themselves, which would be needed to compare peaks properly for ii in range(270, 280): plt.plot(data[ii] - peds[ii], drawstyle="steps", label="pix{}".format(ii)) plt.legend() .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_008.png :alt: raw data exploration :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 234-245 Camera Displays --------------- It’s of course much easier to see the signal if we plot it in 2D with correct pixel positions! note: the instrument data model is not fully implemented, so there is not a good way to load all the camera information (right now it is hacked into the ``inst`` sub-container that is read from the Monte-Carlo file) .. GENERATED FROM PYTHON SOURCE LINES 245-248 .. code-block:: Python :lineno-start: 246 camgeom = source.subarray.tel[24].camera.geometry .. GENERATED FROM PYTHON SOURCE LINES 249-257 .. code-block:: Python :lineno-start: 249 title = "CT24, run {} event {} ped-sub".format(event.index.obs_id, event.index.event_id) disp = CameraDisplay(camgeom, title=title) disp.image = sums - peds disp.cmap = plt.cm.RdBu_r disp.add_colorbar() disp.set_limits_percent(95) # autoscale .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_009.png :alt: CT24, run 1 event 5101 ped-sub :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_009.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 258-266 It looks like a nice signal! We have plotted our pedestal-subtracted trace integral, and see the shower clearly! Let’s look at all telescopes: note we plot here the raw signal, since we have not calculated the pedestals for each) .. GENERATED FROM PYTHON SOURCE LINES 266-280 .. code-block:: Python :lineno-start: 267 for tel in event.r0.tel.keys(): plt.figure() camgeom = source.subarray.tel[tel].camera.geometry title = "CT{}, run {} event {}".format( tel, event.index.obs_id, event.index.event_id ) disp = CameraDisplay(camgeom, title=title) disp.image = event.r0.tel[tel].waveform[0].sum(axis=1) disp.cmap = plt.cm.RdBu_r disp.add_colorbar() disp.set_limits_percent(95) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_010.png :alt: CT26, run 1 event 5101 :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_010.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_011.png :alt: CT27, run 1 event 5101 :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_011.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_012.png :alt: CT29, run 1 event 5101 :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_012.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_013.png :alt: CT47, run 1 event 5101 :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_013.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_014.png :alt: CT53, run 1 event 5101 :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_014.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_015.png :alt: CT59, run 1 event 5101 :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_015.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_016.png :alt: CT69, run 1 event 5101 :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_016.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_017.png :alt: CT73, run 1 event 5101 :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_017.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_018.png :alt: CT121, run 1 event 5101 :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_018.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_019.png :alt: CT122, run 1 event 5101 :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_019.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_020.png :alt: CT124, run 1 event 5101 :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_020.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_021.png :alt: CT148, run 1 event 5101 :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_021.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_022.png :alt: CT162, run 1 event 5101 :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_022.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_023.png :alt: CT166, run 1 event 5101 :srcset: /auto_examples/tutorials/images/sphx_glr_raw_data_exploration_023.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 281-287 some signal processing… ----------------------- Let’s try to detect the peak using the scipy.signal package: https://docs.scipy.org/doc/scipy/reference/signal.html .. GENERATED FROM PYTHON SOURCE LINES 287-305 .. code-block:: Python :lineno-start: 289 pix_ids = np.arange(len(data)) has_signal = sums > 300 widths = np.array( [ 8, ] ) # peak widths to search for (let's fix it at 8 samples, about the width of the peak) peaks = [signal.find_peaks_cwt(trace, widths) for trace in data[has_signal]] for p, s in zip(pix_ids[has_signal], peaks): print("pix{} has peaks at sample {}".format(p, s)) plt.plot(data[p], drawstyle="steps-mid") plt.scatter(np.array(s), data[p, s]) .. GENERATED FROM PYTHON SOURCE LINES 306-309 clearly the signal needs to be filtered first, or an appropriate wavelet used, but the idea is nice .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 7.266 seconds) .. _sphx_glr_download_auto_examples_tutorials_raw_data_exploration.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: raw_data_exploration.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: raw_data_exploration.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: raw_data_exploration.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_