Tool to generate a HDF table of sampling interval coefficients #690

SeiyaNozaki · 2021-04-12T08:53:47Z

This PR is a continuation of #255.
I added a new tool and class to create sampling interval coefficients tables, which is used for TimeSamplingCorrection.

This tool has three stages:

count pulse peak on each capacitor
- we usually use 1 million events, so better to submit a subrun-wise job to save time. It seems some modules have some problems with this method (almost synchronization between test pulse and the first capacitor id), but this trouble appears randomly for all of the modules. So we take 2-3 runs for a single gain to generate sampling interval coefficients.
stack single fits file, convert into sampling interval coefficients, calculate charge resolution with new coefficients for the verification
- stack multiple fits files generated by subrun-wise jobs and convert the stucked peak counts into sampling interval coefficients for each run. Then, charge resolutions are calculated with new coefficients for the verification. If a given module has trouble in a specific run, charge resolution would be not good than expected. So the final coefficient values would be what achieved the best resolutions. The final coefficients and used run for the generation of coefficients on each module are saved for each gain.
merge coefficients of both gains into a single HDF table
- merge output files of the previous stage for each gain. The generated HDF file should be compatible with TimeSamplingCorrection already implemented in lstchain.
  https://github.com/cta-observatory/cta-lstchain/blob/master/lstchain/calib/camera/time_sampling_correction.py#L14

…fficient

maxnoe · 2021-04-12T09:05:05Z

lstchain/calib/camera/sampling_interval_coefficient_calculate.py

+        The SamplingIntervalCalculate class to create a sampling interval coefficient table for LST readout system using chip DRS4.
+    """
+    def __init__(self):
+        self.peak_count = np.zeros([N_PIXELS, N_CAPACITORS_CHANNEL])


zeros creates a float64 array by default, counts should be integers, right? So specify dtype=np.int64 or another appropriate int dtype.

thanks, I set the proper data type (uint16 is enough).

maxnoe · 2021-04-12T09:05:12Z

lstchain/calib/camera/sampling_interval_coefficient_calculate.py

+    """
+    def __init__(self):
+        self.peak_count = np.zeros([N_PIXELS, N_CAPACITORS_CHANNEL])
+        self.fc_count = np.zeros([N_PIXELS, N_CAPACITORS_CHANNEL])


maxnoe · 2021-04-12T09:05:23Z

lstchain/calib/camera/sampling_interval_coefficient_calculate.py

+        self.charge_reso_array_after_corr ={}
+
+        self.charge_reso_final = np.zeros(N_PIXELS)
+        self.used_run = np.zeros(N_PIXELS)


Run is also integer

maxnoe · 2021-04-12T09:07:06Z

lstchain/calib/camera/sampling_interval_coefficient_calculate.py

+            self.fc_count[np.arange(N_PIXELS), fc % N_CAPACITORS_CHANNEL] += 1
+
+
+    def stuck_single_sampling_interval(self, file_list, gain):


What does stuck mean here? Are you sure that's the right word?

sorry, it is a typo... it should be 'stack'

maxnoe · 2021-04-12T09:07:33Z

lstchain/calib/camera/sampling_interval_coefficient_calculate.py

+
+            if gain_in_file == gain:
+                if not run_id in self.peak_count_stuck.keys():
+                    self.peak_count_stuck[run_id] = np.zeros([N_PIXELS, N_CAPACITORS_CHANNEL])


maxnoe · 2021-04-12T09:07:44Z

lstchain/calib/camera/sampling_interval_coefficient_calculate.py

+            if gain_in_file == gain:
+                if not run_id in self.peak_count_stuck.keys():
+                    self.peak_count_stuck[run_id] = np.zeros([N_PIXELS, N_CAPACITORS_CHANNEL])
+                    self.fc_count_stuck[run_id] = np.zeros([N_PIXELS, N_CAPACITORS_CHANNEL])


codecov · 2021-04-12T09:09:11Z

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 85.05%. Comparing base (640eb41) to head (51970ce).
Report is 2679 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #690      +/-   ##
==========================================
+ Coverage   84.35%   85.05%   +0.69%     
==========================================
  Files          61       61              
  Lines        5004     5285     +281     
==========================================
+ Hits         4221     4495     +274     
- Misses        783      790       +7

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

maxnoe · 2021-04-12T09:09:20Z

lstchain/calib/camera/sampling_interval_coefficient_calculate.py

+            self.sampling_interval_coefficient[run_id] = np.zeros([N_PIXELS, N_CAPACITORS_CHANNEL + N_SAMPLES])
+
+            for pixel in range(N_PIXELS):
+                self.sampling_interval_coefficient[run_id][pixel, :N_CAPACITORS_CHANNEL] = \


Please use implicit line continuations inside parentheses instead of these line breaks with \.

Also, it looks like this loop is not necessary but could be avoided by using the axis keyword for the numpy functions.

thanks for the comment. done.
(Concerning the loop, please let me know if there is a more straightforward way...!)

maxnoe · 2021-04-12T13:17:34Z

lstchain/calib/camera/sampling_interval_coefficient_calculate.py

+        r1_sample_end = r0_r1_calibrator.r1_sample_end.tel[tel_id]
+
+        # Check pulse
+        pulse_event_flag = np.sum(np.max(waveform[gain], axis=1) > 100) > 1800


This needs more context and probably the unnamed constants here should be transformed into either global constants or configurable traitlets.

added contexts, and set global constants for them

maxnoe · 2021-04-12T13:17:50Z

lstchain/calib/camera/sampling_interval_coefficient_calculate.py

+
+
+    def set_charge_array(self, gain):
+        self.charge_array_before_corr = np.zeros([N_PIXELS, 60000])


I modified related codes. If max_events is given, this given value will be used to generate a new array.
If not (max_events is None by default), just use a large enough number (55000) to cover all events in a single subrun.
(A single subrun usually contains ~53000 events.)
I'm sure it is not the best way... so I'm very happy if you have any good idea...

You can actually get the number of events by doing len(source.multi_file), I will also add __len__ to the LSTEventSource, but that will require a new release.

A more flexible solution would be to use a python list and append to it, converting to a 2d numpy array at the end.

maxnoe · 2021-04-12T13:17:56Z

lstchain/calib/camera/sampling_interval_coefficient_calculate.py

+        self.charge_reso_array_before_corr = np.zeros(N_PIXELS)
+
+        for run_id in self.peak_count_stuck.keys():
+            self.charge_array_after_corr[run_id] = np.zeros([N_PIXELS, 60000])


same as above

maxnoe · 2021-04-12T13:18:17Z

lstchain/calib/camera/sampling_interval_coefficient_calculate.py

+
+__all__ = ['SamplingIntervalCalculate']
+
+N_PIXELS = 1855


you can import these from ctapipe_io_lst.constants.

maxnoe · 2021-04-12T13:23:22Z

lstchain/tools/lstchain_create_sampling_interval_coefficient_file.py

+        self.log.debug('start peak counting')
+
+        for i, event in enumerate(self.eventsource):
+            if event.index.event_id % 500 == 0:


Use tqdm for a proper progress bar instead.

ok, I didn't understand tqdm before:) Now tqdm is used during even loop.

Then this line is not needed anymore

maxnoe · 2021-04-13T12:47:10Z

lstchain/tools/lstchain_create_sampling_interval_coefficient_file.py

+from lstchain.calib.camera.sampling_interval_coefficient_calculate import SamplingIntervalCalculate
+from lstchain.paths import run_info_from_filename
+
+class SamplingIntervalCoefficientHDFWriter(Tool):


This seems just to be three tools in one that don't really share code and even use different traitlets. I think really doing three smaller tools would be better and less confusing (since the help will only show the options for each step, not for all steps).

I see. Actually, I just wanted to gather all of the procedures in a single file. Is it possible to put three smaller tools in a single file and call one of the tools when executing the command...? or one file per one tool?

You can do multiple Tools per file and have different entrypoints, but I would just make three files

FrancaCassol · 2021-04-27T15:24:50Z

lstchain/tools/lstchain_create_sampling_interval_coefficient_file.py

+    lstchain_create_sampling_interval_coefficient_file --help
+    """
+
+    name = "SamplingIntervalCoefficientHDFWriter"


Hi Seiya,
if I understand correctly you are writing fits files, instead of a HDF5 file.
Would be it possible to write indeed a h5 file? (we decided to prefer this format)

Forget my comment, I see that the final file is in h5 and also that you put together the gains at the ends.
Perhaps a bit misleading the fact that you have single tool for the three steps

…tion to plot results)

SeiyaNozaki · 2021-05-01T08:42:50Z

Thank you for the comments, @maxnoe and @FrancaCassol !
I refactored the codes to make them simple and easy to understand.
Now there are two steps (lstchain_create_peak_count_file_for_sampling_interval_calib and lstchain_create_sampling_interval_coefficient_file): the first one is to count the pulse peak on each capacitor and the second is to stack, convert and verify the sampling interval coefficients.
The summary of the procedure is below:

I also added a function to produce a result pdf file like below:
sampling_interval_coefficient_results_HG.pdf

It would be very appreciated to have a look again.

SeiyaNozaki · 2021-08-04T03:39:06Z

Hi @FrancaCassol @maxnoe,
Sorry for the delay, but came back again.
I'm happy if you can check the updated codes again when you have time!

FrancaCassol · 2021-08-06T13:38:17Z

lstchain/calib/camera/sampling_interval_coefficient_calculate.py

+    def plot_results(self, verify_data_path, output_file):
+
+        plt.rcParams['font.size']=15
+        camera = CameraGeometry.from_name('LSTCam-002')


This is an outdated camera (we have now version 4).
Better to use:

from ctapipe_io_lst import load_camera_geometry camera = load_camera_geometry() camera = camera.transform_to(EngineeringCameraFrame())

FrancaCassol · 2021-08-06T14:02:48Z

lstchain/tools/lstchain_create_sampling_interval_coefficient_file.py

-    input_fits = traits.Path(
-        help="Path to the generated peak count fits files",
+    input_peak_count_fits = traits.Path(
+        help="Path to the generated peak count fits files by lstchain_create_peak_count_file_for_sampling_interval_calib",


Hi @SeiyaNozaki ,

I do not see this tool in your PR, am I wrong?

maxnoe · 2022-01-31T08:56:51Z

lstchain/tools/lstchain_create_sampling_interval_coefficient_file.py

+        self.sampling_interval_coefficient_calculate.convert_to_samp_interval_coefficient()
+        self.sampling_interval_coefficient_calculate.set_charge_array()
+
+        for i, event in enumerate(tqdm(


No need for enumerate, the event already has event.count

SeiyaNozaki added 2 commits April 12, 2021 07:54

add a tool to create sampling interval coefficients

ed75d28

Merge remote-tracking branch 'upstream/master' into samp_interval_coe…

12c3556

…fficient

maxnoe reviewed Apr 12, 2021

View reviewed changes

SeiyaNozaki added 7 commits April 13, 2021 02:31

fix typo

15a78d9

set uint16 as dtype for some parameters

da1ce4a

import ctapipe_io_lst.constants for global parameters

fbf1660

add general constants and explanations

4ae4602

remove line breaks with backslash

b5255df

use tqdm

89ca030

import tqdm

8f23ed7

maxnoe reviewed Apr 13, 2021

View reviewed changes

FrancaCassol reviewed Apr 27, 2021

View reviewed changes

refactoring (separate procedures into two Tools and files, add a func…

51970ce

…tion to plot results)

FrancaCassol reviewed Aug 6, 2021

View reviewed changes

maxnoe reviewed Jan 31, 2022

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Tool to generate a HDF table of sampling interval coefficients #690

Tool to generate a HDF table of sampling interval coefficients #690

SeiyaNozaki commented Apr 12, 2021 •

edited

Loading

maxnoe Apr 12, 2021 •

edited

Loading

SeiyaNozaki Apr 13, 2021

maxnoe Apr 12, 2021

maxnoe Apr 12, 2021

maxnoe Apr 12, 2021

SeiyaNozaki Apr 12, 2021

maxnoe Apr 12, 2021 •

edited

Loading

maxnoe Apr 12, 2021

codecov bot commented Apr 12, 2021 •

edited

Loading

maxnoe Apr 12, 2021

SeiyaNozaki Apr 13, 2021

maxnoe Apr 12, 2021

SeiyaNozaki Apr 13, 2021

maxnoe Apr 12, 2021

SeiyaNozaki Apr 13, 2021

maxnoe Apr 13, 2021

maxnoe Apr 12, 2021

maxnoe Apr 12, 2021

SeiyaNozaki Apr 13, 2021

maxnoe Apr 12, 2021

SeiyaNozaki Apr 13, 2021

maxnoe Jan 31, 2022

maxnoe Apr 13, 2021

SeiyaNozaki Apr 16, 2021

maxnoe Apr 27, 2021

FrancaCassol Apr 27, 2021

FrancaCassol Apr 27, 2021

SeiyaNozaki commented May 1, 2021

SeiyaNozaki commented Aug 4, 2021

FrancaCassol Aug 6, 2021

FrancaCassol Aug 6, 2021

maxnoe Jan 31, 2022 •

edited

Loading

		self.fc_count[np.arange(N_PIXELS), fc % N_CAPACITORS_CHANNEL] += 1


		def stuck_single_sampling_interval(self, file_list, gain):



		def set_charge_array(self, gain):
		self.charge_array_before_corr = np.zeros([N_PIXELS, 60000])


		__all__ = ['SamplingIntervalCalculate']

		N_PIXELS = 1855

Tool to generate a HDF table of sampling interval coefficients #690

Are you sure you want to change the base?

Tool to generate a HDF table of sampling interval coefficients #690

Conversation

SeiyaNozaki commented Apr 12, 2021 • edited Loading

maxnoe Apr 12, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

maxnoe Apr 12, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

codecov bot commented Apr 12, 2021 • edited Loading

Codecov Report

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

SeiyaNozaki commented May 1, 2021

SeiyaNozaki commented Aug 4, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

maxnoe Jan 31, 2022 • edited Loading

Choose a reason for hiding this comment

SeiyaNozaki commented Apr 12, 2021 •

edited

Loading

maxnoe Apr 12, 2021 •

edited

Loading

maxnoe Apr 12, 2021 •

edited

Loading

codecov bot commented Apr 12, 2021 •

edited

Loading

maxnoe Jan 31, 2022 •

edited

Loading