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

muon analysis r0 to dl1 specific function #556

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open

muon analysis r0 to dl1 specific function #556

Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
4375f16
move muon analysis in r0 to dl1 in its own function
vuillaut Nov 5, 2020
b3e19ac
mirror area is already a Quantity
vuillaut Nov 5, 2020
920efe7
docstring
vuillaut Nov 5, 2020
336a6a2
Merge branch 'master' into clean_muon_anaysis_r0_dl1
vuillaut Nov 8, 2020
49a2e1d
Merge branch 'master' into clean_muon_anaysis_r0_dl1
vuillaut Dec 23, 2020
6361784
Merge branch 'master' into clean_muon_anaysis_r0_dl1
rlopezcoto Feb 10, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
185 changes: 106 additions & 79 deletions lstchain/reco/r0_to_dl1.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -238,19 +238,13 @@ def r0_to_dl1(

cal_mc = load_calibrator_from_config(config, subarray)

# minimum number of pe in a pixel to include it
# in calculation of muon ring time (peak sample):
min_pe_for_muon_t_calc = 10.

# Dictionary to store muon ring parameters
muon_parameters = create_muon_table()

if not is_simu:

# TODO : add DRS4 calibration config in config file, read it and pass it here
r0_r1_calibrator = LSTR0Corrections(
pedestal_path=pedestal_path, tel_id=1,
)
r0_r1_calibrator = LSTR0Corrections(pedestal_path=pedestal_path, tel_id=1)

# all this will be cleaned up in a next PR related to the configuration files

Expand Down Expand Up @@ -640,9 +634,6 @@ def r0_to_dl1(
dl1_container.trigger_time = event.r0.tel[telescope_id].trigger_time
dl1_container.trigger_type = event.r0.tel[telescope_id].trigger_type

# FIXME: no need to read telescope characteristics like foclen for every event!
foclen = subarray.tel[telescope_id].optics.equivalent_focal_length
mirror_area = u.Quantity(subarray.tel[telescope_id].optics.mirror_area, u.m ** 2)
dl1_container.prefix = tel.prefix

# extra info for the image table
Expand All @@ -661,75 +652,13 @@ def r0_to_dl1(

# Muon ring analysis, for real data only (MC is done starting from DL1 files)
if not is_simu:
bad_pixels = event.mon.tel[telescope_id].calibration.unusable_pixels[0]
# Set to 0 unreliable pixels:
image = tel.image*(~bad_pixels)

# process only promising events, in terms of # of pixels with large signals:
if tag_pix_thr(image):

# re-calibrate r1 to obtain new dl1, using a more adequate pulse integrator for muon rings
numsamples = event.r1.tel[telescope_id].waveform.shape[2] # not necessarily the same as in r0!
bad_pixels_hg = event.mon.tel[telescope_id].calibration.unusable_pixels[0]
bad_pixels_lg = event.mon.tel[telescope_id].calibration.unusable_pixels[1]
# Now set to 0 all samples in unreliable pixels. Important for global peak
# integrator in case of crazy pixels! TBD: can this be done in a simpler
# way?
bad_waveform = np.array(([np.transpose(np.array(numsamples*[bad_pixels_hg])),
np.transpose(np.array(numsamples*[bad_pixels_lg]))]))

# print('hg bad pixels:',np.where(bad_pixels_hg))
# print('lg bad pixels:',np.where(bad_pixels_lg))

event.r1.tel[telescope_id].waveform *= ~bad_waveform
r1_dl1_calibrator_for_muon_rings(event)

tel = event.dl1.tel[telescope_id]
image = tel.image*(~bad_pixels)

# Check again: with the extractor for muon rings (most likely GlobalPeakWindowSum)
# perhaps the event is no longer promising (e.g. if it has a large time evolution)
if not tag_pix_thr(image):
good_ring = False
else:
# read geometry from event.inst. But not needed for every event. FIXME?
geom = subarray.tel[telescope_id].\
camera.geometry

muonintensityparam, dist_mask, \
ring_size, size_outside_ring, muonringparam, \
good_ring, radial_distribution, \
mean_pixel_charge_around_ring,\
muonpars = \
analyze_muon_event(subarray,
event.index.event_id,
image, geom, foclen,
mirror_area, False, '')
# mirror_area, True, './')
# (test) plot muon rings as png files

# Now we want to obtain the waveform sample (in HG and LG) at which the ring light peaks:
bright_pixels_waveforms = event.r1.tel[telescope_id].waveform[:, image > min_pe_for_muon_t_calc, :]
stacked_waveforms = np.sum(bright_pixels_waveforms, axis=-2)
# stacked waveforms from all bright pixels; shape (ngains, nsamples)
hg_peak_sample = np.argmax(stacked_waveforms, axis=-1)[0]
lg_peak_sample = np.argmax(stacked_waveforms, axis=-1)[1]

if good_ring:
fill_muon_event(-1,
muon_parameters,
good_ring,
event.index.event_id,
dragon_time,
muonintensityparam,
dist_mask,
muonringparam,
radial_distribution,
ring_size,
size_outside_ring,
mean_pixel_charge_around_ring,
muonpars,
hg_peak_sample, lg_peak_sample)
r0_dl1_muon_analysis(event,
telescope_id,
r1_dl1_calibrator_for_muon_rings,
subarray,
muon_parameters,
dragon_time,
)

# writes mc information per telescope, including photo electron image
if is_simu \
Expand Down Expand Up @@ -838,3 +767,101 @@ def rescale_dl1_charge(event, scaling_factor):

for tel_id, tel in event.dl1.tel.items():
tel.image *= scaling_factor


def r0_dl1_muon_analysis(event, telescope_id, r1_dl1_calibrator_for_muon_rings, subarray, muon_parameters, dragon_time):
"""
Muon analysis in the stage r0 to dl1.
Fills the `muon_parameters` table.

Parameters
----------
event: `ctapipe.containers.DataContainer`
telescope_id: int
r1_dl1_calibrator_for_muon_rings: `LSTCameraCalibrator`
subarray: `ctapipe.instrument.subarray.SubarrayDescription`
muon_parameters: dict `lstchain.image.muon.create_muon_table`
dragon_time: `astropy.Quantity`

"""
# minimum number of pe in a pixel to include it
# in calculation of muon ring time (peak sample):
min_pe_for_muon_t_calc = 10.

# FIXME: no need to read telescope characteristics like foclen for every event!
Copy link
Collaborator

Choose a reason for hiding this comment

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

We can remove this FIXME comment, it is a pretty old comment in which I misinterpreted that this would imply some re-reading from disk.

foclen = subarray.tel[telescope_id].optics.equivalent_focal_length
mirror_area = subarray.tel[telescope_id].optics.mirror_area

bad_pixels = event.mon.tel[telescope_id].calibration.unusable_pixels[0]
# Set to 0 unreliable pixels:
image = event.dl1.tel[telescope_id].image*(~bad_pixels)
Copy link
Collaborator

Choose a reason for hiding this comment

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

Is there in ctapipe already some pixel interpolator? Would be interesting to test it. Should be better in case of isolated bad pixels (perhaps a bit dangerous when a whole pixel cluster is unreliable).

Copy link
Member

Choose a reason for hiding this comment

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

I don't think there is, but yes, I think it is necessary.

Copy link
Member

Choose a reason for hiding this comment

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

It should be easy enough to implement approaches like a mean of correct neighbors using the neighbor matrix.


# process only promising events, in terms of # of pixels with large signals:
if tag_pix_thr(image):

# re-calibrate r1 to obtain new dl1, using a more adequate pulse integrator for muon rings
numsamples = event.r1.tel[telescope_id].waveform.shape[2] # not necessarily the same as in r0!
bad_pixels_hg = event.mon.tel[telescope_id].calibration.unusable_pixels[0]
bad_pixels_lg = event.mon.tel[telescope_id].calibration.unusable_pixels[1]
# Now set to 0 all samples in unreliable pixels. Important for global peak
# integrator in case of crazy pixels! TBD: can this be done in a simpler
# way?
bad_waveform = np.array(([np.transpose(np.array(numsamples*[bad_pixels_hg])),
np.transpose(np.array(numsamples*[bad_pixels_lg]))]))

# print('hg bad pixels:',np.where(bad_pixels_hg))
# print('lg bad pixels:',np.where(bad_pixels_lg))

event.r1.tel[telescope_id].waveform *= ~bad_waveform
r1_dl1_calibrator_for_muon_rings(event)

tel = event.dl1.tel[telescope_id]
image = tel.image*(~bad_pixels)

# Check again: with the extractor for muon rings (most likely GlobalPeakWindowSum)
# perhaps the event is no longer promising (e.g. if it has a large time evolution)
if not tag_pix_thr(image):
good_ring = False
else:
# read geometry from event.inst. But not needed for every event. FIXME?
camera_geometry = subarray.tel[telescope_id].camera.geometry

muonintensityparam, dist_mask, \
ring_size, size_outside_ring, muonringparam, \
good_ring, radial_distribution, \
mean_pixel_charge_around_ring, \
muonpars = analyze_muon_event(subarray,
event.index.event_id,
image, camera_geometry,
foclen,
mirror_area,
False,
'',
)
# mirror_area, True, './')
# (test) plot muon rings as png files

# Now we want to obtain the waveform sample (in HG and LG) at which the ring light peaks:
bright_pixels_waveforms = event.r1.tel[telescope_id].waveform[:, image > min_pe_for_muon_t_calc, :]
stacked_waveforms = np.sum(bright_pixels_waveforms, axis=-2)
# stacked waveforms from all bright pixels; shape (ngains, nsamples)
hg_peak_sample = np.argmax(stacked_waveforms, axis=-1)[0]
lg_peak_sample = np.argmax(stacked_waveforms, axis=-1)[1]

if good_ring:
fill_muon_event(-1,
muon_parameters,
good_ring,
event.index.event_id,
dragon_time,
muonintensityparam,
dist_mask,
muonringparam,
radial_distribution,
ring_size,
size_outside_ring,
mean_pixel_charge_around_ring,
muonpars,
hg_peak_sample,
lg_peak_sample,
)