Save pixels charge and peak time arrays in 'DL1' h5 tables

[Elisa-Visentin]

Hi,
Sorry for disturbing you again (and for the question, which is a bit 'crazy').
I need to save the pixel-wise charge and peak time information in DL1 data, preferentially in the same h5 table where Hillas parameters are saved. I managed to extract, for every image, two arrays, one for the charge, one for the peak time (these arrays are not the original ones, too; I had to set at zero some of the values). I used DL1Container to save, for every telescope, the DL1CameraContainer where image and peak_time are saved, but I can't save these arrays in the h5 cells (with HDF5TableWriter). I would like to have an output table where, for every event_id, obs_id and tel_id, both single values (Hillas parameters) and arrays (pixels charge) are saved.
I cannot use cta-process and I can only use ctapipe containers. And I know this is not consistent with the official DL1 data format, but I need to propagate pixel-wise information and save it in a table preserving ctapipe containers.
Is there a way to produce this table?
Thank you

[vuillaut]

Hi @Elisa-Visentin

This is not "crazy" at all and is even already supported.
Why can't you use ctapipe-process? It already writes images and peak charges in the HDF5 files (see --write-images option).

[Elisa-Visentin]

Since which version has ctapipe-process being implemented? I'm using an old ctapipe version (0.12) and I cannot upgrade it now. Then there are other issues with ctapipe-process, which have already been reported... Anyway, I cannot use it 😭

[kosack]

This should be no problem, even in an older version of ctapipe. You simply need to fill a Container with that information and pass it to a TableWriter, which can take a list of containers to write to a single table (all Containers in the list get combined into a single output table). The only problem you might have when writing these to a single table is that the length of each table cell must be the same always, meaning that you can't write e.g. 2000 pixels to one cell, and then a different number of pixels to the next: the TableWriter assumes fixed-length arrays (for reasons of speed optimization). I only mention that because it could be the cause of a problem if you try to write the images from mixed telescope types to the same table.

As it sounds like you did this, can you explain better what you tried? And also explain better what didn't work? It's hard to give a good answer without more information.

[Elisa-Visentin]

The code is something like:

with HDF5TableWriter(output_file, group_name="events", mode="w") as writer:
               
                ....
                
		camera_geom_masked = camera_geoms[tel_id][signal_pixels]
                image_masked = image[signal_pixels]
                peak_time_masked = peak_time[signal_pixels]
                not_signal_pixels=[not item for item in signal_pixels]
                image_all_masked=np.where(not_signal_pixels, 0, image)   #array of charge in every pixel
                time_all=np.where(not_signal_pixels, 0, peak_time) #same with peak time
                
                ....
                
                dl1tel=DL1Container()
                dl1tel.tel[tel_id]=DL1CameraContainer(image=(image_all_masked.astype('float32')),peak_time=(time_all.astype('float32')))
              	
              	...
              	
              	event=EventIndexContainer(
              		obs_id=event.index.obs_id, 
              		event_id=event.index.event_id
                )
                
                writer.write(
                    "image",
                    (event, dl1tel)
            	)

I cut it to remove all the unnecessary steps (cleaning, Hillas parameters...).
This way I would have Hillas and image values in two different tables (a bit annoying but I can merge tables on the event indexes). But in the table I only have two columns (event and obs ID) and no charge and time values. In the log, there was a warning concerning the inability to save the DL1Container in the file.

For each event, and for each telescope, the array length is the same (~8000 pixels) because I set to 0 the ones not surviving cleaning instead of deleting them from the array.

If you need more information, I will try to provide it (sorry, they asked me to do this and I never departed from our standard data format before, so I must disturb you)

[maxnoe]

the TableWriter assumes fixed-length arrays (for reasons of speed optimization).

Not really. It's much more fundamental than just an assumption for speed: the underlying library pytables does not support variable length arrays in structured hdf5 tables.

[maxnoe]

The code is something like:

You call the DL1Container, which is a container for all telescopes dl1tel and set one of its members to a DL1CameraContainer.

What you need to pass to the writer is directly the DL1CameraContainer, you need to write it to a table per telescope (or at least per telescope type) due to the fixed length limitation.

So the correct code would look something like this:

dl1_tel = DL1CameraContainer(image=..., peak_time=...)
writer.write((index, dl1_tel), f"dl1/event/telescope/images/tel_{tel_id:03d}")

This is all handled in recent versions of ctapipe in the DataWriter class.

[Elisa-Visentin]

Thank you, this seems working (at least, now, arrays are in the tables; I must still work on retrieving them). Sorry, I was trying to save them where they could not be saved🤦

BTW, please try to preserve back compatibility while producing new releases: we could already use ctapipe 0.19 but we would have no LST data😅😅

[kosack]

BTW, please try to preserve back compatibility while producing new releases: we could already use ctapipe 0.19 but we would have no LST data😅😅

Yes the purpose of introducing the ctapipe.io.DataWriter (introduced in ctapipe-0.12) and the ctapipe.io.HDF5EventSource was exactly to define a standard data format for DL1 and DL2 information, so that we don't end up with different implementations for each telescope. That way we can transparently handle changes to the data model written to output files. However, LSTChain started a bit too early to benefit from that, but I think should be migrated soon to use the standard ctapipe-process tool which will solve this problem in the future. Note that the DataWriter in the version of ctapipe you have already works and does what you implemented above, so you could simply use it. The format it writes is slightly different from the more recent versions, as the data model has evolved, but it would still make your code easier to maintain.

ctapipe-process defines not only a standard output format (because it uses a DataWriter internally), but also defines a standard analysis and configuration that works for any IACT, with the full mixed-telescope array of CTAO as the main goal of course. This replaces a large amount of code that is currently in lstchain, nectarchain, and protopipe for example, and if used, makes it much easier for us to make changes without breaking things.

[kosack]

Just to re-iterate, both the DataWriter class and the ctapipe-process tool exist in ctapipe-v0.12, so you might save time using them already. At least swapping a hand-rolled HDF5TableWriter (tablewriter(table, containers)) with DataWriter (writer(event)) would simplify your problem down to one line.

see the code in ctapipe/tools/process.py as an example.

[Elisa-Visentin]

Yes, but there are still other issues with ctapipe-process (MAGIC cleaning...), so we need to implement our custom steps.

Then, are there some huge differences between, e.g., Subarray versions (1.0, 2.0)? Keeping as a compatible one even 1.0 would solve half of our issues

[maxnoe]

We cannot, at this point in ctapipe development, provide backwards compatibility guarantees, that would dramatically limit our ability to develop the new functionality we still need or fix bugs in an appropriate way.

This is very clearly put at the top of the README, please understand that.

We also cannot provide a lot of support for relatively old versions of ctapipe, especially since most problems you mentioned are already fixed in newer versions.

There will be more breaking changes ahead, especially since CTAO now has several data model specifiacations released, to which we need to adapt.

[Elisa-Visentin]

We should not use older versions (I agree on this: I spent a lot of nights upgrading the pipeline instead of sleeping) but we cannot use newer versions and data formats. Nice 'impasse'...

[kosack]

Once we have a 1.0 release, though, the idea is to have long-term support of both data format and APIs. We are getting closer to that, but still have a few major changes to implement. So for now, as Max said, we can't really guarantee anything other than starting with the raw data, which should always work, even if newer versions of ctapipe can't read data written by older lstchain versions. Another option would be to implement an EventSource for lstchain-DL1 data, but it's probably not worth the effort as even LST will have a new DL0 format once it is operated by CTAO.

[maxnoe]

To add some more lstchain perspective:

• Update of lstchain to ctapipe 0.19 is underway and is almost finished, however, for now we just updated the ctapipe version, not switched to the ctapipe data format for dl1 and dl2
• There are still many custom steps performed in lstchain, that are not easy to integrate into ctapipe-process, many need updates to ctapipe to do properly
• There are custom fields in the lst data (now as debug in the CTAO R1 data model) , that at least for the commissioning are written and needed for calibration
• I am not sure if the same is true for MAGIC data

But even just extending the DataWriter with a couple of methods to also write the custom data would already help a lot vs. completely starting from scratch.