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.

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

Lint code #74

Merged
merged 13 commits into from
Oct 27, 2023
1 change: 1 addition & 0 deletions .pre-commit-config.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ repos:
hooks:
- id: flake8
args: [--max-line-length=88, "--extend-ignore=E203,E712"]
exclude: ^src/nectarchain/user_scripts/
- repo: https://github.com/mwouts/jupytext
rev: v1.14.6
hooks:
Expand Down
159 changes: 99 additions & 60 deletions examples/waveform_animation_highGain.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,39 +2,43 @@
Script showing how to generate movies with the waveform evolution.

To be used with https://github.com/cta-observatory/ctapipe,
as well as https://github.com/cta-observatory/ctapipe_io_nectarcam, which includes ctapipe.io.nectarcameventsource,
to read NectarCAM data, pending https://github.com/cta-observatory/ctapipe_io_nectarcam/pull/2
as well as https://github.com/cta-observatory/ctapipe_io_nectarcam, which includes
ctapipe.io.nectarcameventsource, to read NectarCAM data, pending
https://github.com/cta-observatory/ctapipe_io_nectarcam/pull/2

protozfitsreader should be installed (see https://github.com/cta-sst-1m/protozfitsreader) because
protozfitsreader should be installed (see
https://github.com/cta-sst-1m/protozfitsreader) because
ctapipe.io.nectarcameventsource depends on it to read zfits files.

"""

import numpy as np
import matplotlib.pylab as plt
from matplotlib.animation import FuncAnimation

import numpy as np
from ctapipe.image import HillasParameterizationError, hillas_parameters, tailcuts_clean
from ctapipe.io import EventSeeker, event_source
from ctapipe.visualization import CameraDisplay
from ctapipe.image import tailcuts_clean, hillas_parameters, HillasParameterizationError
from ctapipe.io import event_source, EventSeeker

from matplotlib.animation import FuncAnimation

ped_path = 'ped/NectarCAM.Run1247.0000.fits.fz'
ped_path = "ped/NectarCAM.Run1247.0000.fits.fz"

num_events = 10

path = 'obs/NectarCAM.Run1250.0000.fits.fz'
path = "obs/NectarCAM.Run1250.0000.fits.fz"
reader = event_source(input_url=path, max_events=num_events)
seeker = EventSeeker(reader)

runid = path.split('NectarCAM.Run')[1].split('.')[0]
runid = path.split("NectarCAM.Run")[1].split(".")[0]

delay = 100 # ms
fps = 1000./delay # frame per second
fps = 1000.0 / delay # frame per second
hi, lo = 0, 1 # channel index in data
tel = 0 # NectarCAM QM telescope ID
adc_to_pe = 58. # theoretical value
winsize_start, winsize_end = 6, 10 # quick rise, long tail: start winsize_start before max, ends winsize_end after max, in ns
adc_to_pe = 58.0 # theoretical value
winsize_start, winsize_end = (
6,
10,
) # quick rise, long tail: start winsize_start before max, ends winsize_end after
# max, in ns
tailcuts = [5, 10]
neighbors = 3 # number of neighboring pixels to include in Hillas cleaning

Expand All @@ -43,7 +47,7 @@
verbose = True
save = False

cmap = 'gnuplot2'
cmap = "gnuplot2"

currentevent = seeker[0]
num_samples = seeker[0].nectarcam.tel[tel].svc.num_samples
Expand All @@ -66,22 +70,22 @@ def get_window(event, gain):
w = event.r0.tel[tel].waveform[gain]
max = np.max(w[w.argmax(axis=0)])
imax = np.where(w.T == max)[0][0]
if imax-winsize_start < 0:
if imax - winsize_start < 0:
istart = 0
iend = winsize_start+winsize_end
elif imax+winsize_end > num_samples-1:
istart = num_samples-1-winsize_start-winsize_end
iend = num_samples-1
iend = winsize_start + winsize_end
elif imax + winsize_end > num_samples - 1:
istart = num_samples - 1 - winsize_start - winsize_end
iend = num_samples - 1
else:
# quick rise, long tail, offset time window
istart = imax-winsize_start
iend = imax+winsize_end
if iend > num_samples-1:
iend = num_samples-1
istart = imax - winsize_start
iend = imax + winsize_end
if iend > num_samples - 1:
iend = num_samples - 1
if istart < 0:
istart = 0
return istart, iend


counter = get_window(currentevent, hi)[0]
hillasdone = False
Expand All @@ -102,7 +106,7 @@ def animation():
disp_hi_raw = CameraDisplay(camgeom, ax=ax_hi_raw, autoupdate=True)
disp_hi_raw.cmap = cmap
disp_hi_raw.add_colorbar(ax=ax_hi_raw)

disp_hi_charge = CameraDisplay(camgeom, ax=ax_hi_charge, autoupdate=True)
disp_hi_charge.cmap = cmap
disp_hi_charge.add_colorbar(ax=ax_hi_charge)
Expand All @@ -111,61 +115,82 @@ def animation():
disp_hi_wf.cmap = cmap
disp_hi_wf.add_colorbar(ax=ax_hi_wf)

ax_hi_plot_wf.set_title('High gain, wave form (ADC)')
ax_hi_plot_wf.set_xlabel('Time (ns)')
ax_hi_plot_wf.set_ylabel('ADC summed over all pixels')
ax_hi_plot_wf.set_title("High gain, wave form (ADC)")
ax_hi_plot_wf.set_xlabel("Time (ns)")
ax_hi_plot_wf.set_ylabel("ADC summed over all pixels")

def update(frames):
global counter, currentevent, hillasdone, sw, ped_hi
event = currentevent

hiistart, hiiend = get_window(event, hi)

if counter >= hiiend or event.r0.tel[tel].trigger_type != trigtype:
event = next(iter(seeker)) # get next event
currentevent = event
disp_hi_charge.clear_overlays()
hiistart, hiiend = get_window(event, hi)
hillasdone = False # reset Hillas reco flag
counter = hiistart

w = event.r0.tel[tel].waveform

ax_hi_raw.set_title("High gain, raw data (ADC), event {}".format(event.r0.event_id))

ax_hi_raw.set_title(
"High gain, raw data (ADC), event {}".format(event.r0.event_id)
)

hiw2 = w[hi].T[hiistart:hiiend] # select time window
wf_hi_max = np.amax(hiw2)

if verbose:
print('INFO High gain: event id {}, counter {}, max(waveform)={}, window {}-{} ns'.format(
event.r0.event_id, counter, w[hi].max(), hiistart, hiiend))

ax_hi_charge.set_title("High gain, charge (PE), event {}, window {}-{} ns".format(
event.r0.event_id, hiistart, hiiend))
print(
"INFO High gain: event id {}, counter {}, max(waveform)={}, window "
"{}-{} ns".format(
event.r0.event_id, counter, w[hi].max(), hiistart, hiiend
)
)

ax_hi_charge.set_title(
"High gain, charge (PE), event {}, window {}-{} ns".format(
event.r0.event_id, hiistart, hiiend
)
)

image_hi_raw = hiw2.sum(axis=0)

if perform_calib:
# Very rough calibration
image_hi_charge = ((w[hi]-peds_hi).T[hiistart:hiiend].T/adc_to_pe).sum(axis=1)
image_hi_charge = ((w[hi] - peds_hi).T[hiistart:hiiend].T / adc_to_pe).sum(
axis=1
)

if clean and not hillasdone:
# Cleaning
cleanmask_hi = tailcuts_clean(camgeom, image_hi_charge,
picture_thresh=tailcuts[1],
boundary_thresh=tailcuts[0],
min_number_picture_neighbors=neighbors)
cleanmask_hi = tailcuts_clean(
camgeom,
image_hi_charge,
picture_thresh=tailcuts[1],
boundary_thresh=tailcuts[0],
min_number_picture_neighbors=neighbors,
)
image_hi_charge[cleanmask_hi == 0] = 0

# Hillas reco
try:
hillas_param_hi = hillas_parameters(camgeom, image_hi_charge)
disp_hi_charge.overlay_moments(hillas_param_hi, with_label=False,
color='red', alpha=0.7,
linewidth=2, linestyle='dashed')
disp_hi_charge.highlight_pixels(cleanmask_hi, color='white', alpha=0.3, linewidth=2)

sw.append(hillas_param_hi.width.value/hillas_param_hi.intensity)
disp_hi_charge.overlay_moments(
hillas_param_hi,
with_label=False,
color="red",
alpha=0.7,
linewidth=2,
linestyle="dashed",
)
disp_hi_charge.highlight_pixels(
cleanmask_hi, color="white", alpha=0.3, linewidth=2
)

sw.append(hillas_param_hi.width.value / hillas_param_hi.intensity)
except HillasParameterizationError:
disp_hi_charge.clear_overlays()
disp_hi_charge.axes.figure.canvas.draw()
Expand All @@ -174,7 +199,7 @@ def update(frames):

charge_hi = image_hi_charge.sum()
if verbose:
print(' charge hi = {} pe'.format(charge_hi))
print(" charge hi = {} pe".format(charge_hi))

disp_hi_raw.image = image_hi_raw
disp_hi_raw.set_limits_percent(95)
Expand All @@ -184,25 +209,39 @@ def update(frames):
disp_hi_charge.set_limits_percent(95)
disp_hi_charge.axes.figure.canvas.draw()

disp_hi_wf.image = hiw2[counter-hiistart]
disp_hi_wf.image = hiw2[counter - hiistart]
disp_hi_wf.set_limits_minmax(0, wf_hi_max)
disp_hi_wf.axes.figure.canvas.draw()

ax_hi_wf.set_title("High gain, wave form (ADC), event {}, time {} ns".format(event.r0.event_id, counter))
ax_hi_wf.set_title(
"High gain, wave form (ADC), event {}, time {} ns".format(
event.r0.event_id, counter
)
)
ax_hi_plot_wf.plot(w[hi].sum(axis=0)[0:counter])
ax_hi_plot_wf.figure.canvas.draw()

counter += 1 # beurk...
return [ax_hi_wf, ax_hi_plot_wf, ax_hi_raw, ax_hi_charge, ]
return [
ax_hi_wf,
ax_hi_plot_wf,
ax_hi_raw,
ax_hi_charge,
]

try:
frames = num_events*(winsize_start+winsize_end)
except:
frames = num_events * (winsize_start + winsize_end)
except Exception:
frames = None
anim = FuncAnimation(fig, update, repeat=False,
interval=delay, frames=frames, blit=(not save))
anim = FuncAnimation(
fig, update, repeat=False, interval=delay, frames=frames, blit=(not save)
)
if save:
anim.save(filename=path.replace('.fits.fz', '_highGain.mp4'), fps=fps, extra_args=['-vcodec', 'libx264'])
anim.save(
filename=path.replace(".fits.fz", "_highGain.mp4"),
fps=fps,
extra_args=["-vcodec", "libx264"],
)
plt.show()


Expand All @@ -211,5 +250,5 @@ def main():
animation()


if __name__ == '__main__':
if __name__ == "__main__":
main()
24 changes: 13 additions & 11 deletions notebooks/Access NectarCAM data using DIRAC API.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,11 @@
# extension: .py
# format_name: percent
# format_version: '1.3'
# jupytext_version: 1.14.5
# jupytext_version: 1.14.6
# kernelspec:
# display_name: Python (nectar-dev)
# display_name: nectarchain-ctapipe0.19
# language: python
# name: nectar-dev
# name: nectarchain-ctapipe0.19
# ---

# %% [markdown]
Expand All @@ -18,37 +18,39 @@
# In this short notebook, we will see how to access data stored on the grid using the
# DIRAC API, without the hassle of first manually downloading them locally.
#
# In order to achieve this, you will obviously need a `conda` environment in which all
# relevant code is installed, such as `ctapipe`, `nectarchain`, but also `CTADIRAC`
# itself. Please refer to the `nectarchain` installation procedure at
# <https://github.com/cta-observatory/nectarchain> and follow the instructions to enable
# DIRAC support.
# In order to achieve this, you will obviously need a `conda` environment in which
# all relevant code is installed, such as `ctapipe`, `nectarchain`, but also
# `CTADIRAC` itself. Please refer to the `nectarchain` installation procedure at
# <https://github.com/cta-observatory/nectarchain> and follow the instructions to
# enable DIRAC support.
#
# You will also need to have an active proxy for EGI, initialized e.g. with:
#
# ```
# dirac-proxy-init -U -M -g cta_nectarcam
# dirac-proxy-init -M -g cta_nectarcam
# ```
#
# You can also check whether you currently have an active proxy with the command
# `dirac-proxy-info`.

# %%
# %matplotlib inline
import os
from glob import glob

from ctapipe.coordinates import EngineeringCameraFrame
from ctapipe.instrument import CameraGeometry
from ctapipe.io import EventSeeker, EventSource
from ctapipe.visualization import CameraDisplay

# %%
# %matplotlib inline
from DIRAC.Interfaces.API.Dirac import Dirac

# %%
dirac = Dirac()

# %%
# ?dirac.getFile
# dirac.getFile?

# %%
lfns = [
Expand Down
8 changes: 5 additions & 3 deletions notebooks/Read_NectarCAM_events.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,9 @@
# name: python3
# ---

# %%
# %matplotlib inline

import time

# %%
Expand Down Expand Up @@ -39,7 +42,6 @@

# %%
evt = seeker.get_event_index(25)

adcsum = evt.r0.tel[0].waveform[0].sum(axis=1)

camera = CameraGeometry.from_name("NectarCam-003")
Expand Down Expand Up @@ -128,8 +130,8 @@
disp.cmap = cmap
disp.add_colorbar()

# Comment: if this cell says that "charges" is not defined it's because the event type
# is not 1. Re-run it.
# Comment: if this cell says that "charges" is not defined it's because the event
# type is not 1. Re-run it.

# %% [markdown]
# ## Hillas cleaning
Expand Down
3 changes: 0 additions & 3 deletions src/nectarchain/data/container/__init__.py
Original file line number Diff line number Diff line change
@@ -1,3 +0,0 @@
from .chargesContainer import *
from .core import *
from .waveformsContainer import *
Loading
Loading