Page for users' guide on pedestal estimation tool (#159)

* first stab at doc for pedestal estimation tool

* links and fixes

* additions/corrections suggested by JPL

Author: tibaldo

docs/user-guide/howto-pedestal.rst

@@ -0,0 +1,62 @@
+.. _pedestal:
+
+Pedestal estimation tool
+-------------------------------
+The pedestal estimation tool processes pedestal calibration data
+(dedicated pedestal runs or runs containing interleaved pedestal
+events) and returns estimates of the pedestal in ADC counts as a function of gain channel/pixel/sample within the readout window.
+
+Run the tool locally
+=========================
+To use the pedestal estimation tool, first activate the ``nectarchain`` environment:
+
+.. code-block:: console
+
+    $ mamba activate nectarchain
+
+
+The user script `nectarchain/user_scripts/ltibaldo/example_pedestal.py` showcases the usage of the tool.
+
+The input data are identified by run number. See :ref:`env-vars` to
+set up the ``$NECTARCAMDATA`` environment variable. Events in the
+selected runs are
+processed in slices with a fixed number of events set by the
+``events_per_slice`` parameter (see `~nectarchain.makers.core.EventsLoopNectarCAMCalibrationTool`).
+
+The pedestal
+estimation tool inherits the configurable parameters of the
+`~nectarchain.makers.component.PedestalComponent.PedestalEstimationComponent`.
+The data can be filtered based on time using the ``ucts_tmin`` and
+``ucts_tmax`` parameters and to eliminate outlier waveforms using the ``filter_method`` parameter. Two different methods to exclude outlier
+waveforms are implemented:
+
+* ``WaveformsStdFilter`` discards waveforms with a standard deviation
+  exceeding the threshold value defined by the parameter
+  ``wfs_std_threshold``;
+  
+* ``ChargeDistributionFilter`` discards waveforms with a total charge integrated over the entire readout window in the tails of the charge distribution, either below ``charge_sigma_low_thr`` or above ``charge_sigma_high_thr``.
+
+
+To run the example script:
+
+.. code-block:: console
+
+    $ python -i example_pedestal.py
+
+Inspect the results
+=========================
+The results are stored in a
+`~nectarchain.data.container.pedestalContainer.NectarCAMPedestalContainer`. The
+results include information on pixels that were flagged as having
+an abnormal behavior during the computation of the pedestals. The
+flags are defined in in
+`~nectarchain.data.container.pedestalContainer.PedestalFlagBits`. The
+results are accessible on the fly if the tool is run interactively (as in the example above) and stored in a `.h5` file.
+
+The user script `nectarchain/user_scripts/ltibaldo/show_pedestal_output.py` provides an example of how to access the results from disk and produce some plots:
+
+.. code-block:: console
+
+    $ python -i plot_pedestal_output.py
+
+

docs/user-guide/index.rst

@@ -10,4 +10,5 @@ User Guide
    env-vars
    howto-dirac
    dqm
+   howto-pedestal
    troubleshooting

src/nectarchain/user_scripts/ltibaldo/example_pedestal.py

@@ -8,14 +8,12 @@
 log = logging.getLogger(__name__)
 log.handlers = logging.getLogger("__main__").handlers
 
-from nectarchain.makers.calibration import (
-    PedestalNectarCAMCalibrationTool,
-)
+from nectarchain.makers.calibration import PedestalNectarCAMCalibrationTool
 
 run_number = 3938
-max_events= 2999
+max_events = 2999
 events_per_slice = 300
-outfile = "/Users/ltibaldo/tmp/test_pedestal/pedestal_{}.h5".format(run_number)
+outfile = os.environ["NECTARCAMDATA"] + "/runs/pedestal_{}.h5".format(run_number)
 
 tool = PedestalNectarCAMCalibrationTool(
     progress_bar=True,
@@ -24,13 +22,13 @@
     events_per_slice=events_per_slice,
     log_level=20,
     output_path=outfile,
-    overwrite = True,
-    filter_method = 'WaveformsStdFilter',
-    wfs_std_threshold = 4.,
+    overwrite=True,
+    filter_method="WaveformsStdFilter",
+    wfs_std_threshold=4.0,
 )
 
 tool.initialize()
 tool.setup()
 
 tool.start()
-output = tool.finish(return_output_component=True)
+output = tool.finish(return_output_component=True)

src/nectarchain/user_scripts/ltibaldo/plot_pedestal_output.py

@@ -1,13 +1,15 @@
+import os
+
+import matplotlib.pyplot as plt
 import numpy as np
 import tables
-import matplotlib.pyplot as plt
 
-filename = '/Users/ltibaldo/tmp/test_pedestal/pedestal_3938.h5'
+filename = os.environ["NECTARCAMDATA"] + "/runs/pedestal_3938.h5"
 pixel_id = 132
 sample = 13
 nsamples = 60
 
-container_name = 'NectarCAMPedestalContainer'
+container_name = "NectarCAMPedestalContainer"
 
 h5file = tables.open_file(filename)
 
@@ -20,41 +22,50 @@
 # fill results for plotting
 i = 0
 for result in h5file.root.__members__:
-    table = h5file.root[result]['NectarCAMPedestalContainer'][0]
-    wf = table['pedestal_mean_hg'][table['pixels_id'] == pixel_id][0]
-    std = table['pedestal_std_hg'][table['pixels_id'] == pixel_id][0]
-    if result == 'data_combined':
+    table = h5file.root[result]["NectarCAMPedestalContainer"][0]
+    wf = table["pedestal_mean_hg"][table["pixels_id"] == pixel_id][0]
+    std = table["pedestal_std_hg"][table["pixels_id"] == pixel_id][0]
+    if result == "data_combined":
         pedestal_combined = wf
         pedestal_std = std
     else:
         pedestals[i] = wf
         pedestals_std[i] = std
-        tmin = np.append(tmin,table['ucts_timestamp_min'])
-        tmax = np.append(tmax, table['ucts_timestamp_max'])
+        tmin = np.append(tmin, table["ucts_timestamp_min"])
+        tmax = np.append(tmax, table["ucts_timestamp_max"])
         i += 1
 
 tmean = 0.5 * (tmin + tmax)
 
 fig1 = plt.figure()
 ax1 = plt.axes()
 ax1.set_title(f"Pixel {pixel_id}")
-ax1.set_xlabel('sample (ns)')
-ax1.set_ylabel('pedestal (ADC counts)')
+ax1.set_xlabel("sample (ns)")
+ax1.set_ylabel("pedestal (ADC counts)")
 for s in range(len(pedestals)):
-    ax1.plot(pedestals[s],color='0.5',alpha=0.2)
-ax1.plot(pedestal_combined,color='r')
+    ax1.plot(pedestals[s], color="0.5", alpha=0.2)
+ax1.plot(pedestal_combined, color="r")
 
 fig2 = plt.figure()
 ax2 = plt.axes()
 ax2.set_title(f"Pixel {pixel_id} Sample {sample}")
-ax2.set_xlabel('UCTS timestamp')
-ax2.set_ylabel('pedestal (ADC counts)')
-ax2.errorbar(tmean,pedestals[:,sample],
-             xerr=[tmean-tmin,tmax-tmean],
-             yerr=[pedestals_std[:,sample]],
-             fmt = 'o', color='k',capsize=0.)
-ax2.axhspan(pedestal_combined[sample]-pedestal_std[sample],
-            pedestal_combined[sample]+pedestal_std[sample],
-            facecolor='r', alpha=0.2, zorder=0)
-
-plt.show()
+ax2.set_xlabel("UCTS timestamp")
+ax2.set_ylabel("pedestal (ADC counts)")
+ax2.errorbar(
+    tmean,
+    pedestals[:, sample],
+    xerr=[tmean - tmin, tmax - tmean],
+    yerr=[pedestals_std[:, sample]],
+    fmt="o",
+    color="k",
+    capsize=0.0,
+)
+ax2.axhspan(
+    pedestal_combined[sample] - pedestal_std[sample],
+    pedestal_combined[sample] + pedestal_std[sample],
+    facecolor="r",
+    alpha=0.2,
+    zorder=0,
+)
+
+plt.show()