Merge pull request #46 from lucabaldini/more_digi

More digi

docs/api.rst

@@ -27,6 +27,7 @@ List of modules
    mc
    modeling
    pprint
+   readout
    recon
    rng
    roi

docs/digi.rst

@@ -1,20 +1,16 @@
 :mod:`hexsample.digi` --- Digitization facilities
 =================================================
 
-This module contains all the digitization facilities, including the event structure
-for a digitized events (:class:`DigiEvent <hexsample.digi.DigiEvent>`) and a generic
-event readout structure (:class:`HexagonalReadout <hexsample.digi.HexagonalReadout>`).
+This module contains...
 
 Digi event structure
 --------------------
 
 .. literalinclude:: ../hexsample/digi.py
-   :pyobject: DigiEvent
+   :pyobject: DigiEventRectangular
    :end-before: __post_init__
 
 
-Digitization process
---------------------
 
 
 

docs/readout.rst

@@ -0,0 +1,15 @@
+:mod:`hexsample.readout` --- Readout facilities
+===============================================
+
+This module contains...
+
+
+Digitization process
+--------------------
+
+
+
+Module documentation
+--------------------
+
+.. automodule:: hexsample.readout

docs/release.rst

@@ -4,6 +4,8 @@ Release notes
 =============
 
 
+* Merging https://github.com/lucabaldini/hexsample/pull/46
+* Implemented circular and sparse readout modes.
 * Minor tweaks to the docs to make clear that all the lengths are in cm.
 * Minor formatting changes.
 

hexsample/app.py

@@ -24,6 +24,7 @@
 
 from hexsample import __pkgname__, __version__, __tagdate__, __url__
 from hexsample.hexagon import HexagonalLayout
+from hexsample.readout import HexagonalReadoutMode
 
 
 START_MESSAGE = f"""
@@ -154,6 +155,11 @@ def add_readout_options(self) -> None:
             help='number of rows in the readout chip')
         group.add_argument('--pitch', type=float, default=0.005,
             help='pitch of the readout chip in cm')
+        modes = [item.value for item in HexagonalReadoutMode]
+        group.add_argument('--readoutmode', type=str, choices=modes, default='RECTANGULAR',
+            help='readout mode')
+        group.add_argument('--padding', type=int, nargs=4, default=(2, 2, 2, 2),
+            help='padding on the four sides of the ROT')
         group.add_argument('--noise', type=float, default=20.,
             help='equivalent noise charge rms in electrons')
         group.add_argument('--gain', type=float, default=1.,
@@ -164,8 +170,6 @@ def add_readout_options(self) -> None:
             help='trigger threshold in electron equivalent')
         group.add_argument('--zsupthreshold', type=int, default=0,
             help='zero-suppression threshold in ADC counts')
-        group.add_argument('--padding', type=int, nargs=4, default=(2, 2, 2, 2),
-            help='padding on the four sides of the ROT')
 
     def add_sensor_options(self) -> None:
         """Add an option group for the sensor properties.

hexsample/bin/hxdisplay.py

@@ -24,7 +24,7 @@
 
 from hexsample import logger
 from hexsample.app import ArgumentParser
-from hexsample.digi import HexagonalReadout
+from hexsample.readout import HexagonalReadoutRectangular
 from hexsample.display import HexagonalGridDisplay
 from hexsample.fileio import DigiInputFile
 from hexsample.hexagon import HexagonalLayout
@@ -48,7 +48,7 @@ def hxdisplay(**kwargs):
     header = input_file.header
     args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\
         header['pitch'], header['noise'], header['gain']
-    readout = HexagonalReadout(*args)
+    readout = HexagonalReadoutRectangular(*args)
     logger.info(f'Readout chip: {readout}')
     display = HexagonalGridDisplay(readout)
     for event in input_file:

hexsample/bin/hxplots.py

@@ -0,0 +1,86 @@
+#!/usr/bin/env python
+#
+# Copyright (C) 2023 luca.baldini@pi.infn.it
+#
+# For the license terms see the file LICENSE, distributed along with this
+# software.
+#
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, write to the Free Software Foundation Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+"""Event file viewer comparing reconstructed quantities with MC truth.
+"""
+from ast import literal_eval
+
+import numpy as np
+
+from hexsample.app import ArgumentParser
+from hexsample.fileio import ReconInputFile
+from hexsample.hist import Histogram1d, Histogram2d
+from hexsample.plot import plt, setup_gca
+from hexsample.analysis import create_histogram
+
+
+__description__ = \
+"""Simple viewer for comparing reconstructed energy and position with the MC
+true values.
+"""
+
+# Parser object.
+HXVIEW_ARGPARSER = ArgumentParser(description=__description__)
+HXVIEW_ARGPARSER.add_infile()
+
+def hxview(**kwargs):
+    """View the file content.
+    Shows histograms of energy and cluster_size of recon events vs their MC truth.
+    """
+    input_file = ReconInputFile(kwargs['infile'])
+    # Plotting the reconstructed energy and the true energy
+    histo = create_histogram(input_file, 'energy', mc = False)
+    mc_histo = create_histogram(input_file, 'energy', mc = True)
+    plt.figure('Photons energy')
+    histo.plot(label='Reconstructed')
+    mc_histo.plot(label='MonteCarlo')
+    plt.xlabel('Energy [eV]')
+    plt.legend()
+
+    # Plotting the reconstructed x and y position and the true position.
+    plt.figure('Reconstructed photons position')
+    binning = np.linspace(-5. * 0.1, 5. * 0.1, 100)
+    x = input_file.column('posx')
+    y = input_file.column('posy')
+    histo = Histogram2d(binning, binning).fill(x, y)
+    histo.plot()
+    setup_gca(xlabel='x [cm]', ylabel='y [cm]')
+    plt.figure('True photons position')
+    x_mc = input_file.mc_column('absx')
+    y_mc = input_file.mc_column('absy')
+    histo_mc = Histogram2d(binning, binning).fill(x_mc, y_mc)
+    histo_mc.plot()
+    setup_gca(xlabel='x [cm]', ylabel='y [cm]')
+    #Closing the file and showing the figures.
+    plt.figure('x-direction resolution')
+    binning = np.linspace((x-x_mc).min(), (x-x_mc).max(), 100)
+    histx = Histogram1d(binning, xlabel=r'$x - x_{MC}$ [cm]').fill(x-x_mc)
+    histx.plot()
+    plt.figure('y-direction resolution')
+    binning = np.linspace((y-y_mc).min(), (y-y_mc).max(), 100)
+    histy = Histogram1d(binning, xlabel=r'$y - y_{MC}$ [cm]').fill(y-y_mc)
+    histy.plot()
+
+    input_file.close()
+    plt.show()
+
+if __name__ == '__main__':
+    hxview(**vars(HXVIEW_ARGPARSER.parse_args()))

hexsample/bin/hxrecon.py

@@ -27,8 +27,8 @@
 from hexsample import logger
 from hexsample.app import ArgumentParser, check_required_args
 from hexsample.clustering import ClusteringNN
-from hexsample.digi import HexagonalReadout
-from hexsample.fileio import DigiInputFile, ReconOutputFile
+from hexsample.readout import HexagonalReadoutMode, HexagonalReadoutSparse, HexagonalReadoutRectangular, HexagonalReadoutCircular
+from hexsample.fileio import DigiInputFileBase, DigiInputFileSparse, DigiInputFileRectangular, DigiInputFileCircular, DigiInputFile, ReconOutputFile, peek_readout_type
 from hexsample.hexagon import HexagonalLayout
 from hexsample.recon import ReconEvent
 
@@ -52,21 +52,55 @@ def hxrecon(**kwargs):
     input_file_path = str(kwargs['infile'])
     if not input_file_path.endswith('.h5'):
         raise RuntimeError('Input file {input_file_path} does not look like a HDF5 file')
-    input_file = DigiInputFile(input_file_path)
-    header = input_file.header
-    args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\
-        header['pitch'], header['noise'], header['gain']
-    readout = HexagonalReadout(*args)
-    logger.info(f'Readout chip: {readout}')
+
+    # It is necessary to extract the reaodut type because every readout type
+    # corresponds to a different DigiEvent type.
+    # If there is no info about the readout, we try to reconstruct with a Rectangular
+    # readout mode, that is the mode for all the old reconstruction before the 
+    # distinction between different readouts was implemented.
+    try:
+        readout_mode = peek_readout_type(input_file_path)
+        # Now we can construct a set of if/else.
+        if readout_mode is HexagonalReadoutMode.SPARSE:
+            input_file = DigiInputFileSparse(input_file_path)
+            header = input_file.header
+            args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\
+                header['pitch'], header['noise'], header['gain']
+            readout = HexagonalReadoutSparse(*args)
+            logger.info(f'Readout chip: {readout}')
+        elif readout_mode is HexagonalReadoutMode.RECTANGULAR:
+            input_file = DigiInputFile(input_file_path)
+            header = input_file.header
+            args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\
+                header['pitch'], header['noise'], header['gain']
+            readout = HexagonalReadoutRectangular(*args)
+            logger.info(f'Readout chip: {readout}')
+        elif readout_mode is HexagonalReadoutMode.CIRCULAR:
+            input_file = DigiInputFileCircular(input_file_path)
+            header = input_file.header
+            args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\
+                header['pitch'], header['noise'], header['gain']
+            readout = HexagonalReadoutCircular(*args)
+            logger.info(f'Readout chip: {readout}')
+    except RuntimeError:
+        input_file = DigiInputFile(input_file_path)
+        header = input_file.header
+        args = HexagonalLayout(header['layout']), header['numcolumns'], header['numrows'],\
+            header['pitch'], header['noise'], header['gain']
+        #readout = HexagonalReadoutCircular(*args)
+        readout = HexagonalReadoutRectangular(*args)
+        logger.info(f'Readout chip: {readout}')
+    # When the readout tipology is determined, the event is clustered ...   
     clustering = ClusteringNN(readout, kwargs['zsupthreshold'], kwargs['nneighbors'])
     suffix = kwargs['suffix']
     output_file_path = input_file_path.replace('.h5', f'_{suffix}.h5')
+    # ... and saved into an output file.
     output_file = ReconOutputFile(output_file_path)
     output_file.update_header(**kwargs)
     output_file.update_digi_header(**input_file.header)
     for i, event in tqdm(enumerate(input_file)):
         cluster = clustering.run(event)
-        args = event.trigger_id, event.timestamp(), event.livetime, event.roi.size, cluster
+        args = event.trigger_id, event.timestamp(), event.livetime, cluster
         recon_event = ReconEvent(*args)
         mc_event = input_file.mc_event(i)
         output_file.add_row(recon_event, mc_event)

hexsample/bin/hxsim.py

@@ -29,8 +29,9 @@
 from hexsample import rng
 from hexsample import HEXSAMPLE_DATA
 from hexsample.app import ArgumentParser
-from hexsample.digi import HexagonalReadout
-from hexsample.fileio import DigiOutputFile
+from hexsample.readout import HexagonalReadoutMode, readout_chip
+from hexsample.fileio import DigiDescriptionSparse, DigiDescriptionRectangular,\
+    DigiDescriptionCircular, digioutput_class
 from hexsample.hexagon import HexagonalLayout
 from hexsample.mc import PhotonList
 from hexsample.roi import Padding
@@ -63,15 +64,24 @@ def hxsim(**kwargs):
     material = Material(kwargs['actmedium'], kwargs['fano'])
     sensor = Sensor(material, kwargs['thickness'], kwargs['transdiffsigma'])
     photon_list = PhotonList(source, sensor, kwargs['numevents'])
+    readout_mode = HexagonalReadoutMode(kwargs['readoutmode'])
+    # Is there any nicer way to do this? See https://github.com/lucabaldini/hexsample/issues/51
+    if readout_mode is HexagonalReadoutMode.SPARSE:
+        readout_args = kwargs['trgthreshold'], kwargs['zsupthreshold'], kwargs['offset']
+    elif readout_mode is HexagonalReadoutMode.RECTANGULAR:
+        padding = Padding(*kwargs['padding'])
+        readout_args = kwargs['trgthreshold'], padding, kwargs['zsupthreshold'], kwargs['offset']
+    elif readout_mode is HexagonalReadoutMode.CIRCULAR:
+        readout_args = kwargs['trgthreshold'], kwargs['zsupthreshold'], kwargs['offset']
+    else:
+        raise RuntimeError
     args = HexagonalLayout(kwargs['layout']), kwargs['numcolumns'], kwargs['numrows'],\
         kwargs['pitch'], kwargs['noise'], kwargs['gain']
-    readout = HexagonalReadout(*args)
+    readout = readout_chip(readout_mode, *args)
     logger.info(f'Readout chip: {readout}')
     output_file_path = kwargs.get('outfile')
-    output_file = DigiOutputFile(output_file_path)
+    output_file = digioutput_class(readout_mode)(output_file_path)
     output_file.update_header(**kwargs)
-    padding = Padding(*kwargs['padding'])
-    readout_args = kwargs['trgthreshold'], padding, kwargs['zsupthreshold'], kwargs['offset']
     logger.info('Starting the event loop...')
     for mc_event in tqdm(photon_list):
         x, y = mc_event.propagate(sensor.trans_diffusion_sigma)

hexsample/bin/hxview.py

@@ -40,8 +40,8 @@
 HXVIEW_ARGPARSER.add_infile()
 HXVIEW_ARGPARSER.add_argument("attribute", type=str, help='Attribute to be viewed.\
                                To be taken from\ the following list:\
-                               trigger_id (int),  timestamp (float),  livetime (int)  \
-                               roi_size (int), energy (float), position (Tuple[float,float])\
+                               trigger_id (int),  timestamp (float),  livetime (int)  \
+                               roi_size (int), energy (float), position (Tuple[float,float])\
                                cluster_size (int), roi_size (int)')
 HXVIEW_ARGPARSER.add_argument("mc_table", type=str, help='Tells if the quantities are in mc table\
                               accepts True or False.')
@@ -62,4 +62,4 @@ def hxview(**kwargs):
     plt.show()
 
 if __name__ == '__main__':
-    hxview(**vars(HXVIEW_ARGPARSER.parse_args()))
+    hxview(**vars(HXVIEW_ARGPARSER.parse_args()))