Merge pull request #209 from pawel21/improve_create_pedestal_script

Improve create pedestal script

lstchain/calib/camera/drs4.py

@@ -1,6 +1,9 @@
 import numpy as np
 from numba import jit, prange
 
+from ctapipe.core import Component
+from ctapipe.core.traits import Int
+
 
 __all__ = ['DragonPedestal']
 
@@ -13,23 +16,32 @@
 roi_size = 40
 high_gain = 0
 low_gain = 1
+n_module_in_camera = 265
 n_gain = 2
 n_channel = 7
 
 
-class DragonPedestal:
+class DragonPedestal(Component):
     """
         The DragonPedestal class to create pedestal
         for LST readout system using chip DRS4.
     """
 
-    def __init__(self, tel_id, n_module):
+    r0_sample_start = Int(default_value=11,
+                          help='Start sample for waveform'
+                          ).tag(config=True)
+
+    def __init__(self, tel_id, n_module, **kwargs):
+        super().__init__(**kwargs)
         self.tel_id = tel_id
-        self.n_module = n_module
-        self.n_pixels = n_module*n_channel # Each module has 7 channels (pixels)
+        self.n_module = n_module # This is number of module read from data
+
+        # Readout system of LST has 265 modules.
+        # Each module has 7 channels (pixels)
+        self.n_pixels = n_module_in_camera*n_channel
         self.meanped = np.zeros((n_gain, self.n_pixels, size4drs))
         self.numped = np.zeros((n_gain, self.n_pixels, size4drs))
-        self.first_cap_array = np.zeros((n_module, n_gain, n_channel))
+        self.first_cap_array = np.zeros((self.n_module, n_gain, n_channel))
 
     def fill_pedestal_event(self, event):
         expected_pixel_id = event.lst.tel[self.tel_id].svc.pixel_ids
@@ -42,26 +54,29 @@ def fill_pedestal_event(self, event):
                                       self.first_cap_array,
                                       self.meanped,
                                       self.numped,
-                                      self.n_module)
+                                      self.n_module,
+                                      self.r0_sample_start)
 
     @staticmethod
     @jit(parallel=True)
     def _fill_pedestal_event_jit(waveform, expected_pixel_id, first_cap_array,
-                                 meanped, numped, n_module):
+                                 meanped, numped, n_module, start_sample_r0):
         for nr_module in prange(0, n_module):
             first_cap = first_cap_array[nr_module, :, :]
             for gain in prange(0, n_gain):
                 for pix in prange(0, n_channel):
                     fc = first_cap[gain, pix]
                     pixel = expected_pixel_id[nr_module * 7 + pix]
-
-                    posads0 = int((2 + fc) % size4drs)
-                    if posads0 + 40 < size4drs:
-                        meanped[gain, pixel, posads0:(posads0 + 36)] += waveform[gain, pixel, 2:38]
-                        numped[gain, pixel, posads0:(posads0 + 36)] += 1
-
+                    posads0 = int((start_sample_r0+fc)%size4drs)
+                    if posads0 + roi_size < size4drs:
+                        # the first 9 samples have occasionally increased signal due to Tsutomu pattern,
+                        # hence we skip them. Start sample might be set as script argument. Default = 11.
+                        meanped[gain, pixel, posads0:((posads0-start_sample_r0) + roi_size-2)] += waveform[gain, pixel, start_sample_r0:roi_size-2]
+                        numped[gain, pixel, posads0:((posads0-start_sample_r0) + roi_size-2)] += 1
                     else:
-                        for k in prange(2, roi_size - 2):
+                        for k in prange(start_sample_r0, roi_size - 2):
+                            # the first 9 samples have occasionally increased signal due to Tsutomu pattern,
+                            # hence we skip them. Start sample might be set as script argument. Default = 11.
                             posads = int((k + fc) % size4drs)
                             val = waveform[gain, pixel, k]
                             meanped[gain, pixel, posads] += val

scripts/lstchain_data_create_pedestal_file.py

@@ -15,24 +15,34 @@
 To run script in console:
 python lstchain_data_create_pedestal_file.py --input_file LST-1.1.Run00097.0000.fits.fz --output_file pedestal.fits 
 --max_events 9000
+To set start sample in waveform --start_r0_waveform i (default i = 11)
 not to use deltaT correction add --deltaT False
 '''
 
 parser = argparse.ArgumentParser()
 
 # Required arguments
-parser.add_argument("--input_file", help="Path to fitz.fz file to create pedestal file.",
-                    type=str, default="")
+parser.add_argument("--input_file",
+                    help="Path to fitz.fz file to create pedestal file.",
+                    type=str)
 
-parser.add_argument("--output_file", help="Path where script create pedestal file",
+parser.add_argument("--output_file",
+                    help="Path where script create pedestal file",
                     type=str)
 
 # Optional argument
-parser.add_argument("--max_events", help="Maximum numbers of events to read."
-                                         "Default = 5000",
-                    type=int, default=5000)
-
-parser.add_argument('--deltaT', '-s', type=lambda x: bool(strtobool(x)),
+parser.add_argument("--max_events",
+                    help="Maximum numbers of events to read. Default = 8000",
+                    type=int,
+                    default=8000)
+
+parser.add_argument("--start_r0_waveform",
+                    help="Start sample for waveform. Default = 11",
+                    type=int,
+                    default=11)
+
+parser.add_argument('--deltaT', '-s',
+                    type=lambda x: bool(strtobool(x)),
                     help='Boolean. True for use deltaT correction'
                     'Default=True, use False otherwise',
                     default=True)
@@ -47,17 +57,17 @@
 
     seeker = EventSeeker(reader)
     ev = seeker[0]
-    telid = ev.r0.tels_with_data[0]
-    n_modules = ev.lst.tel[telid].svc.num_modules
+    tel_id = ev.r0.tels_with_data[0]
+    n_modules = ev.lst.tel[tel_id].svc.num_modules
 
     config = Config({
         "LSTR0Corrections": {
-            "tel_id": telid
+            "tel_id": tel_id
         }
     })
     lst_r0 = LSTR0Corrections(config=config)
 
-    pedestal = DragonPedestal(tel_id=telid, n_module=n_modules)
+    pedestal = DragonPedestal(tel_id=tel_id, n_module=n_modules, r0_sample_start=args.start_r0_waveform)
 
     if args.deltaT:
         print("DeltaT correction active")
@@ -67,7 +77,6 @@
                 pedestal.fill_pedestal_event(event)
                 if i%500 == 0:
                     print("i = {}, ev id = {}".format(i, event.r0.event_id))
-
     else:
         print("DeltaT correction no active")
         for i, event in enumerate(reader):
@@ -78,7 +87,7 @@
     # Finalize pedestal and write to fits file
     pedestal.finalize_pedestal()
 
-    primaryhdu = fits.PrimaryHDU(ev.lst.tel[telid].svc.pixel_ids)
+    primaryhdu = fits.PrimaryHDU(ev.lst.tel[tel_id].svc.pixel_ids)
     secondhdu = fits.ImageHDU(np.int16(pedestal.meanped))
 
     hdulist = fits.HDUList([primaryhdu, secondhdu])