Dqm update - using standard ctapipe tools (#154)

* modified charge intergration using ctapipe charge extraction methods

* adjusting defaults

* fixing camera plotting in charge_integration

* switching from LocalPeakWindowSum to GlobalPeakWindowSum as default

* Simplify existing unit tests

---------

Co-authored-by: jlenain <[email protected]>

src/nectarchain/dqm/camera_monitoring.py

@@ -39,7 +39,7 @@ def __init__(self, gaink):
         self.ChargeInt_Figures_Dict = {}
         self.ChargeInt_Figures_Names_Dict = {}
 
-    def ConfigureForRun(self, path, Pix, Samp, Reader1):
+    def ConfigureForRun(self, path, Pix, Samp, Reader1, **kwargs):
         # define number of pixels and samples
         self.Pix = Pix
         self.Samp = Samp

src/nectarchain/dqm/charge_integration.py

@@ -1,11 +1,21 @@
 import ctapipe.instrument.camera.readout
 import numpy as np
 from ctapipe.coordinates import EngineeringCameraFrame
-from ctapipe.image import LocalPeakWindowSum
+from ctapipe.image.extractor import FixedWindowSum  # noqa: F401
+from ctapipe.image.extractor import FullWaveformSum  # noqa: F401
+from ctapipe.image.extractor import LocalPeakWindowSum  # noqa: F401
+from ctapipe.image.extractor import NeighborPeakWindowSum  # noqa: F401
+from ctapipe.image.extractor import SlidingWindowMaxSum  # noqa: F401
+from ctapipe.image.extractor import TwoPassWindowSum  # noqa: F401
+from ctapipe.image.extractor import GlobalPeakWindowSum
 from ctapipe.visualization import CameraDisplay
+from ctapipe_io_nectarcam import constants
 from matplotlib import pyplot as plt
 from traitlets.config.loader import Config
 
+from ..makers.component import ChargesComponent
+from ..makers.component.core import ArrayDataComponent
+from ..makers.extractor.utils import CtapipeExtractor
 from .dqm_summary_processor import DQMSummary
 
 __all__ = ["ChargeIntegrationHighLowGain"]
@@ -51,7 +61,7 @@ def __init__(self, gaink):
         self.ChargeInt_Figures_Dict = {}
         self.ChargeInt_Figures_Names_Dict = {}
 
-    def ConfigureForRun(self, path, Pix, Samp, Reader1):
+    def ConfigureForRun(self, path, Pix, Samp, Reader1, **charges_kwargs):
         # define number of pixels and samples
         self.Pix = Pix
         self.Samp = Samp
@@ -72,44 +82,60 @@ def ConfigureForRun(self, path, Pix, Samp, Reader1):
         ].camera.readout = ctapipe.instrument.camera.readout.CameraReadout.from_name(
             "NectarCam"
         )
-        config = Config({"LocalPeakWindowSum": {"window_shift": 4, "window_width": 12}})
-
-        self.integrator = LocalPeakWindowSum(subarray, config=config)
+        if charges_kwargs:
+            extractor_kwargs = (
+                ChargesComponent._get_extractor_kwargs_from_method_and_kwargs(
+                    method=charges_kwargs["method"],
+                    kwargs=charges_kwargs["extractor_kwargs"],
+                )
+            )
+            self.integrator = eval(charges_kwargs["method"])(
+                subarray, **extractor_kwargs
+            )
+        else:
+            config = Config(
+                {"GlobalPeakWindowSum": {"window_shift": 4, "window_width": 12}}
+            )
+            self.integrator = GlobalPeakWindowSum(subarray, config=config)
 
     def ProcessEvent(self, evt, noped):
-        self.pixelBAD = evt.mon.tel[0].pixel_status.hardware_failing_pixels
         pixel = evt.nectarcam.tel[0].svc.pixel_ids
-        if len(pixel) < self.Pix:
-            pixel_masked_shutter = list(
-                np.arange(0, self.Pix - len(pixel), 1, dtype=int)
-            )
-            pixel = list(pixel)
-            pixels = np.concatenate([pixel_masked_shutter, pixel])
-        else:
-            pixels = pixel
+        self.pixelBADplot = evt.mon.tel[0].pixel_status.hardware_failing_pixels
+        pixels = pixel
+        self.pixels = pixels
 
-        waveform = evt.r0.tel[0].waveform[self.k]
+        (
+            broken_pixels_hg,
+            broken_pixels_lg,
+        ) = ArrayDataComponent._compute_broken_pixels_event(evt, pixels)
+
+        if self.k == 0:
+            self.pixelBAD = broken_pixels_hg
+            channel = constants.HIGH_GAIN
+        if self.k == 1:
+            self.pixelBAD = broken_pixels_lg
+            channel = constants.LOW_GAIN
 
+        waveform = evt.r0.tel[0].waveform[self.k]
+        waveform = waveform[pixels]
         ped = np.mean(waveform[:, 20])
 
         if noped:
             w_noped = waveform - ped
-            output = self.integrator(
-                w_noped, 0, np.zeros(self.Pix, dtype=int), self.pixelBAD
+            output = CtapipeExtractor.get_image_peak_time(
+                self.integrator(w_noped, 0, channel, self.pixelBAD)
             )
-            image = output.image
-            peakpos = output.peak_time
-            image = image[pixels]
-            peakpos = peakpos[pixels]
+
+            image = output[0]
+            peakpos = output[1]
 
         else:
-            output = self.integrator(
-                waveform, 0, np.zeros(self.Pix, dtype=int), self.pixelBAD
+            output = CtapipeExtractor.get_image_peak_time(
+                self.integrator(waveform, 0, channel, self.pixelBAD)
             )
-            image = output.image
-            peakpos = output.peak_time
-            image = image[pixels]
-            peakpos = peakpos[pixels]
+
+            image = output[0]
+            peakpos = output[1]
 
         if evt.trigger.event_type.value == 32:  # count peds
             self.counter_ped += 1
@@ -274,9 +300,9 @@ def PlotResults(self, name, FigPath):
         # Charge integration MEAN plot
         if self.counter_evt > 0:
             fig1, disp = plt.subplots()
-            disp = CameraDisplay(self.camera[~self.pixelBAD[0]])
+            disp = CameraDisplay(self.camera[~self.pixelBADplot[0]])
             # disp = CameraDisplay(self.subarray.tels[0].camera)
-            disp.image = self.image_all_average[~self.pixelBAD[0]]
+            disp.image = self.image_all_average
             disp.cmap = plt.cm.coolwarm
             disp.axes.text(
                 2,
@@ -302,8 +328,8 @@ def PlotResults(self, name, FigPath):
 
         if self.counter_ped > 0:
             fig2, disp = plt.subplots()
-            disp = CameraDisplay(self.camera[~self.pixelBAD[0]])
-            disp.image = self.image_ped_average[~self.pixelBAD[0]]
+            disp = CameraDisplay(self.camera[~self.pixelBADplot[0]])
+            disp.image = self.image_ped_average
             disp.cmap = plt.cm.coolwarm
             disp.axes.text(
                 2,
@@ -330,8 +356,8 @@ def PlotResults(self, name, FigPath):
         # Charge integration MEDIAN plot
         if self.counter_evt > 0:
             fig3, disp = plt.subplots()
-            disp = CameraDisplay(self.camera[~self.pixelBAD[0]])
-            disp.image = self.image_all_median[~self.pixelBAD[0]]
+            disp = CameraDisplay(self.camera[~self.pixelBADplot[0]])
+            disp.image = self.image_all_median
             disp.cmap = plt.cm.coolwarm
             disp.axes.text(
                 2,
@@ -357,8 +383,8 @@ def PlotResults(self, name, FigPath):
 
         if self.counter_ped > 0:
             fig4, disp = plt.subplots()
-            disp = CameraDisplay(self.camera[~self.pixelBAD[0]])
-            disp.image = self.image_ped_median[~self.pixelBAD[0]]
+            disp = CameraDisplay(self.camera[~self.pixelBADplot[0]])
+            disp.image = self.image_ped_median
             disp.cmap = plt.cm.coolwarm
             disp.axes.text(
                 2,
@@ -385,8 +411,8 @@ def PlotResults(self, name, FigPath):
         # Charge integration STD plot
         if self.counter_evt > 0:
             fig5, disp = plt.subplots()
-            disp = CameraDisplay(self.camera[~self.pixelBAD[0]])
-            disp.image = self.image_all_std[~self.pixelBAD[0]]
+            disp = CameraDisplay(self.camera[~self.pixelBADplot[0]])
+            disp.image = self.image_all_std
             disp.cmap = plt.cm.coolwarm
             disp.axes.text(
                 2,
@@ -412,8 +438,8 @@ def PlotResults(self, name, FigPath):
 
         if self.counter_ped > 0:
             fig6, disp = plt.subplots()
-            disp = CameraDisplay(self.camera[~self.pixelBAD[0]])
-            disp.image = self.image_ped_std[~self.pixelBAD[0]]
+            disp = CameraDisplay(self.camera[~self.pixelBADplot[0]])
+            disp.image = self.image_ped_std
             disp.cmap = plt.cm.coolwarm
             disp.axes.text(
                 2,
@@ -440,8 +466,8 @@ def PlotResults(self, name, FigPath):
         # Charge integration RMS plot
         if self.counter_evt > 0:
             fig7, disp = plt.subplots()
-            disp = CameraDisplay(self.camera[~self.pixelBAD[0]])
-            disp.image = self.image_all_rms[~self.pixelBAD[0]]
+            disp = CameraDisplay(self.camera[~self.pixelBADplot[0]])
+            disp.image = self.image_all_rms
             disp.cmap = plt.cm.coolwarm
             disp.axes.text(
                 2,
@@ -467,8 +493,8 @@ def PlotResults(self, name, FigPath):
 
         if self.counter_ped > 0:
             fig8, disp = plt.subplots()
-            disp = CameraDisplay(self.camera[~self.pixelBAD[0]])
-            disp.image = self.image_ped_rms[~self.pixelBAD[0]]
+            disp = CameraDisplay(self.camera[~self.pixelBADplot[0]])
+            disp.image = self.image_ped_rms
             disp.cmap = plt.cm.coolwarm
             disp.axes.text(
                 2,
@@ -495,7 +521,7 @@ def PlotResults(self, name, FigPath):
         # Charge integration SPECTRUM
         if self.counter_evt > 0:
             fig9, disp = plt.subplots()
-            for i in range(self.Pix):
+            for i in range(len(self.pixels)):
                 plt.hist(
                     self.image_all[:, i],
                     100,
@@ -532,7 +558,7 @@ def PlotResults(self, name, FigPath):
 
         if self.counter_ped > 0:
             fig10, disp = plt.subplots()
-            for i in range(self.Pix):
+            for i in range(len(self.pixels)):
                 plt.hist(
                     self.image_ped[:, i],
                     100,

src/nectarchain/dqm/mean_camera_display.py

@@ -31,7 +31,7 @@ def __init__(self, gaink):
         self.MeanCameraDisplay_Figures_Dict = {}
         self.MeanCameraDisplay_Figures_Names_Dict = {}
 
-    def ConfigureForRun(self, path, Pix, Samp, Reader1):
+    def ConfigureForRun(self, path, Pix, Samp, Reader1, **kwargs):
         # define number of pixels and samples
         self.Pix = Pix
         self.Samp = Samp

src/nectarchain/dqm/mean_waveforms.py

@@ -23,7 +23,7 @@ def __init__(self, gaink):
         self.MeanWaveForms_Figures_Dict = {}
         self.MeanWaveForms_Figures_Names_Dict = {}
 
-    def ConfigureForRun(self, path, Pix, Samp, Reader1):
+    def ConfigureForRun(self, path, Pix, Samp, Reader1, **kwargs):
         # define number of pixels and samples
         self.Pix = Pix
         self.Samp = Samp

src/nectarchain/dqm/pixel_participation.py

@@ -25,7 +25,7 @@ def __init__(self, gaink):
         self.PixelParticipation_Figures_Dict = {}
         self.PixelParticipation_Figures_Names_Dict = {}
 
-    def ConfigureForRun(self, path, Pix, Samp, Reader1):
+    def ConfigureForRun(self, path, Pix, Samp, Reader1, **kwargs):
         # define number of pixels and samples
         self.Pix = Pix
         self.Samp = Samp

src/nectarchain/dqm/pixel_timeline.py

@@ -24,7 +24,7 @@ def __init__(self, gaink):
         self.PixelTimeline_Figures_Dict = {}
         self.PixelTimeline_Figures_Names_Dict = {}
 
-    def ConfigureForRun(self, path, Pix, Samp, Reader1):
+    def ConfigureForRun(self, path, Pix, Samp, Reader1, **kwargs):
         # define number of pixels and samples
         self.Pix = Pix
         self.Samp = Samp

src/nectarchain/dqm/start_dqm.py

@@ -1,4 +1,5 @@
 import argparse
+import json
 import os
 import sys
 import time
@@ -18,6 +19,7 @@
 from nectarchain.dqm.pixel_participation import PixelParticipationHighLowGain
 from nectarchain.dqm.pixel_timeline import PixelTimelineHighLowGain
 from nectarchain.dqm.trigger_statistics import TriggerStatistics
+from nectarchain.makers import ChargesNectarCAMCalibrationTool
 
 
 def main():
@@ -40,6 +42,28 @@ def main():
         action="store_true",
         help="Enables pedestal subtraction in charge integration",
     )
+    # extractor arguments
+    parser.add_argument(
+        "--method",
+        choices=[
+            "FullWaveformSum",
+            "FixedWindowSum",
+            "GlobalPeakWindowSum",
+            "LocalPeakWindowSum",
+            "SlidingWindowMaxSum",
+            "TwoPassWindowSum",
+        ],
+        default="GlobalPeakWindowSum",
+        help="charge extractor method",
+        type=str,
+    )
+    parser.add_argument(
+        "--extractor_kwargs",
+        default='{"window_shift": 4, "window_width": 16}',
+        help="charge extractor kwargs",
+        type=json.loads,
+    )
+
     parser.add_argument(
         "-r",
         "--runnb",
@@ -94,9 +118,21 @@ def main():
     # Defining and printing the options
     PlotFig = args.plot
     noped = args.noped
+    method = args.method
+    extractor_kwargs = args.extractor_kwargs
 
     print("Plot:", PlotFig)
     print("Noped:", noped)
+    print("method:", method)
+    print("extractor_kwargs:", extractor_kwargs)
+
+    kwargs = {"method": method, "extractor_kwargs": extractor_kwargs}
+    charges_kwargs = {}
+    tool = ChargesNectarCAMCalibrationTool()
+    for key in tool.traits().keys():
+        if key in kwargs.keys():
+            charges_kwargs[key] = kwargs[key]
+    print(charges_kwargs)
 
     def GetName(RunFile):
         name = RunFile.split("/")[-1]
@@ -189,7 +225,7 @@ def CreateFigFolder(name, type):
         break
 
     for p in processors:
-        p.ConfigureForRun(path, Pix, Samp, reader1)
+        p.ConfigureForRun(path, Pix, Samp, reader1, **charges_kwargs)
 
     for evt in tqdm(
         reader, total=args.max_events if args.max_events else len(reader), unit="ev"

src/nectarchain/dqm/tests/test_camera_monitoring.py

@@ -9,15 +9,9 @@
 
 
 class TestCameraMonitoring:
-    run_number = 3798
-    max_events = 1
-
     def test_camera_monitoring(self):
-        # run_number = 3938
         path = get_dataset_path("NectarCAM.Run3938.30events.fits.fz")
 
-        config = None
-
         config = Config(
             dict(
                 NectarCAMEventSource=dict(
@@ -29,18 +23,17 @@ def test_camera_monitoring(self):
                 )
             )
         )
-        print(path)
 
         reader1 = EventSource(input_url=path, config=config, max_events=1)
 
         Pix, Samp = CameraMonitoring(HIGH_GAIN).DefineForRun(reader1)
 
-        CameraMonitoring(HIGH_GAIN).ConfigureForRun(path, Pix, Samp, reader1)
-
+        sql_file_date = None
         for evt in tqdm(reader1, total=1):
             run_start1 = evt.nectarcam.tel[0].svc.date
-            SqlFileDate = astropytime.Time(run_start1, format="unix").iso.split(" ")[0]
-        # print("SqlFileDate", SqlFileDate)
-        # CameraMonitoring(HIGH_GAIN).FinishRun()
+            sql_file_date = astropytime.Time(run_start1, format="unix").iso.split(" ")[
+                0
+            ]
+
         assert Pix + Samp == 1915
-        assert SqlFileDate == "2023-01-23"
+        assert sql_file_date == "2023-01-23"