Questions arising during attempt to implement a ReducedWindow and optimize event selection and get IRFs

[mdpunch]

As part of my effort to see if we can mitigate the effect of extreme GRBs, I am trying to see the effect of having a reduced readout window, to reduce the data volume needing to be transferred.

For now, I am just looking at CTAO-N, i.e. MST-NectarCAMs and LSTs. The required readout windows are 60ns and 30ns respectively (see here for justification: "A-PERF-0935 Cherenkov Image Information" https://jama.cta-observatory.org/perspective.req#/items/28805?projectId=11

Reducing the window reduces the impact distance, in effect, so reducing the collection area. But for extreme GRBs we don't care much, so I want to look at reduced windows of 15ns and 20ns for MST-NectarCAMs and LSTs respectively.

First question: how to implement a ReducedWindow:

In discussion with @kosack , he said:

I think the easiest way would be to implement a new DataVolumeReducer subclass that does it. It would just take the new window size as a configuration TelescopeParameter (and anything else you need). Then you should just be able to enable it by setting the config parameters:

CameraCalibrator:
  data_volume_reducer_type: ReadoutWindowReducer

ReadoutWindowReducer:
  max_window_size: 10

The easiest thing would be to add it to the end of ctapipe/src/ctapipe/image/reducer.py which is where the other reducers are.

... if your new subclass is in the code, it should automatically become available to ctapipe-process. If you type "ctapipe-process --help-all", and look in the section called CameraCalibrator. For example see here:

 390   │ --CameraCalibrator.data_volume_reducer_type=<ComponentName>
 391   │     The name of a DataVolumeReducer subclass. Possible values:
 392   │     ['NullDataVolumeReducer', 'TailCutsDataVolumeReducer']
 393   │     Default: 'NullDataVolumeReducer'

If you created a new subclass and put it in the same file as the other reducers (or in a separate file, but one that is imported by the ctapipe.calib module), it should appear there as an option. So thenyou can either set:

--CameraCalibrator.data_volume_reducer_type=YourClassName

or else make a config file with the appropriate option, which would have something like;

CameraCalibrator:
    data_volume_reducer_type: YourClassName

MP (@mdpunch): In the end I just did this the dumb way, defining my own ReadoutWindowReducer(event,subarray) and calling this in a loop prior to the calibrator/image_processor/shower_processor, and writing the DL1 and DL2 to a file. (My main doubt is how the calibrator works on the reduced window... as it's a black box for now ... until I look in the code ; maybe I should apply the calibrator prior to the ReducedWindow until I can solve the problem of calibration on reduced window data?).

Second question: how to apply the new tools for the ML training / Cut optimization / IRF production

Karl gave some links for tutorials for the other tools:

There's even a tutorial now, though it's somewhat generic:

• https://ctapipe.readthedocs.io/en/latest/user-guide/tools/dl2_guide.html
• https://ctapipe.readthedocs.io/en/latest/user-guide/tools/irf_guide.html

My questions on this:

• Anyway, that gives me the gamma/proton/electron files needed for the next step, with and without the reduced window. Following the tutorial, I trained the models, which appeared to work. What's missing for me at this step is some way to visualise the results, for example plots to check for overtraining, plots of the "gammaness" of the test/train samples of each particle, plots of the energy resolution, etc. Is there a notebook somewhere which has examples of this?

• So in the previous step, from six input merged files (3 gamma, 2 proton, 1 electron), I got one file each of {gamma/proton/electron}_final.dl2.h5, but for next step the procedure asks for two of each in input, for the optimization and IRF files respectively. So does that mean I have to use twelve files in input altogether to get two of each now? (since my assumption is that the files which are input here must have the ML applied).

• This leads to my problem now, that there seems to be no ML information in the files (no RandomForest at least). See error message below.

2025-05-02 12:28:52,883 ERROR [ctapipe.ctapipe-optimize-event-selection] (tool.run): Caught unexpected exception: Error evaluating expression '('valid classifier', 'RandomForestClassifier_is_valid')': name 'RandomForestClassifier_is_valid' is not defined
Traceback (most recent call last):
  File "/home/punch/miniconda3/envs/ctapipe_0.24/lib/python3.13/site-packages/ctapipe/core/tool.py", line 445, in run
    self.start()
    ~~~~~~~~~~^^
  File "/home/punch/miniconda3/envs/ctapipe_0.24/lib/python3.13/site-packages/ctapipe/tools/optimize_event_selection.py", line 150, in start
    events, count, meta = loader.load_preselected_events(
                          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^
        self.chunk_size, self.obs_time
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    )
    ^
  File "/home/punch/miniconda3/envs/ctapipe_0.24/lib/python3.13/site-packages/ctapipe/irf/preprocessing.py", line 227, in load_preselected_events
    selected = events[self.epp.quality_query.get_table_mask(events)]
                      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^
  File "/home/punch/miniconda3/envs/ctapipe_0.24/lib/python3.13/site-packages/ctapipe/core/qualityquery.py", line 138, in get_table_mask
    for i, res in enumerate(self.engine(table), start=1):
                  ~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/punch/miniconda3/envs/ctapipe_0.24/lib/python3.13/site-packages/ctapipe/core/expression_engine.py", line 44, in __call__
    raise ExpressionError(
        f"Error evaluating expression '{expression}': {err}"
    ) from None
ctapipe.core.expression_engine.ExpressionError: Error evaluating expression '('valid classifier', 'RandomForestClassifier_is_valid')': name 'RandomForestClassifier_is_valid' is not defined

[maxnoe]

MP (@mdpunch): In the end I just did this the dumb way, defining my own ReadoutWindowReducer(event,subarray) and calling this in a loop prior to the calibrator/image_processor/shower_processor, and writing the DL1 and DL2 to a file. (My main doubt is how the calibrator works on the reduced window... as it's a black box for now ... until I look in the code ; maybe I should apply the calibrator prior to the ReducedWindow until I can solve the problem of calibration on reduced window data?).

I don't think there is a specific issue with the calibrator regarding reduced windows. It would depend on the image extractor you use, but I don't think you use the FullWindowExtractor, right? All peak finding extractors should use a much smaller integration window than what you reduce the waveforms to. I think LST and NectarCam use LocalPeakWindowSum with a width of 8 slices in their analysis. The default in ctapipe is NeighborPeakWindowSum, also with 8 slices window size.

In the current design, you cannot apply the calibrator before the window reducer, as the order of operations is:

R1 to DL0 (this is where data volume reduction is applied, so where your waveform reducer will be run)
DL0 to DL1 (image extraction is run, this is also where calibration coefficients are applied to the resulting charges).

However, the bigger issue with the data volume reducer is that the API limits them to what is currently the assumption of DVR in CTAO: that DVR selects pixels from the waveform data.

You cannot, in the current API, implement a data volume reducer that removes samples from the waveform.

What you could do however, is do it at the ImageExtractor level.

[mdpunch]

MP (@mdpunch): In the end I just did this the dumb way, defining my own ReadoutWindowReducer(event,subarray) and calling this in a loop prior to the calibrator/image_processor/shower_processor, and writing the DL1 and DL2 to a file. (My main doubt is how the calibrator works on the reduced window... as it's a black box for now ... until I look in the code ; maybe I should apply the calibrator prior to the ReducedWindow until I can solve the problem of calibration on reduced window data?).

I don't think there is a specific issue with the calibrator regarding reduced windows. It would depend on the image extractor you use, but I don't think you use the FullWindowExtractor, right?
I just use the supplied yaml files, which indeed have ["type", "*", "NeighborPeakWindowSum"] (nothing special for NectarCAM/LST.

All peak finding extractors should use a much smaller integration window than what you reduce the waveforms to. I think LST and NectarCam use LocalPeakWindowSum with a width of 8 slices in their analysis. The default in ctapipe is NeighborPeakWindowSum, also with 8 slices window size.

This is interesting to know, thanks. For real data, I guess we will use the injected pedestals for the pedestal calibration part (though actually I think we don't do this yet!), but of course for the MCs we don't have these (though maybe it might be interesting to include these in a later prod? ... especially for the high NSB cases...).
I was more worrying that some kind of "image elimination with stand-off" might be needed, to allow the remaining pixels to be used for calibration, which would be quite complicated to implement correctly.
If the remaining slices - after reducing the window to 15ns and removing 8 slices around the peak - are okay for the calibration, that's fine (though it seems a little short...).

In the current design, you cannot apply the calibrator before the window reducer, as the order of operations is:

R1 to DL0 (this is where data volume reduction is applied, so where your waveform reducer will be run) DL0 to DL1 (image extraction is run, this is also where calibration coefficients are applied to the resulting charges).
This is where I put the Reducer routine "by hand".

However, the bigger issue with the data volume reducer is that the API limits them to what is currently the assumption of DVR in CTAO: that DVR selects pixels from the waveform data.

You cannot, in the current API, implement a data volume reducer that removes samples from the waveform.

What you could do however, is do it at the ImageExtractor level.
Of course, if a "GRB mode" gets implemented, it would be the cameras which reduce their windows, so nothing will need to be removed!
But for the simulations, I have to reduce the the samples by hand. Maybe the API could be changed to allow a bit more flexibility?

[maxnoe]

What's missing for me at this step is some way to visualise the results, for example plots to check for overtraining, plots of the "gammaness" of the test/train samples of each particle, plots of the energy resolution, etc. Is there a notebook somewhere which has examples of this?

I have some of this here, from the time the tools were developed: https://github.com/maxnoe/test_ctapipe_ml/blob/main/performance.ipynb

It should also be part of the benchmarking tools that the group of @kosack developed, though I haven't had a detailed look at that yet.

[mdpunch]

I'm running this in 0.24, according to the version when I import from ctapipe.version import version

version 
'0.24.0'

But for now, I can get around this by {not looping over but} summing the n_showers given by read_simulation_configuration()
(possibly having to multiply by shower_reuse)

[mdpunch]

Which works... (though I'll have to get some files like cta_requirements_North-50h.dat, magic_sensitivity_2014.ecsv etc.)

Would you like to look at the outputs?

[maxnoe]

(possibly having to multiply by shower_reuse)

Yes, if you take the n_showers from the simulation config, you have to multiply by shower_reuse.

However, I think we should try to figure out why you don't have the histograms.

[maxnoe]

Would you like to look at the outputs?

Sure, if you want, we can jump into a quick meeting

[mdpunch]

I'm online on Teams.

[maxnoe]

So in the previous step, from six input merged files (3 gamma, 2 proton, 1 electron), I got one file each of {gamma/proton/electron}_final.dl2.h5, but for next step the procedure asks for two of each in input, for the optimization and IRF files respectively. So does that mean I have to use twelve files in input altogether to get two of each now? (since my assumption is that the files which are input here must have the ML applied).

Weren't you the one to remark that we shouldn't optimize the cuts on the same events as we compute the sensitivity? 😉

Anyways, no, you only need one more of each, not twice the amount. You can also in principle, use the same files for each step, but as you correctly remarked, there is a danger of overfitting the cuts on the sensitivity.

So, if you want to use separate datasets for optimizing the cuts from computing sensitivity and IRFs, I think you will want:

• 4 gamma diffuse datasets (train energy regressor, train particle classifier, optimize cuts, sensitivity & IRFs)
• 3 proton datasets (train particle classifier, optimize cuts, sensitivity & IRFs)
• 2 electron datasets (optimize cuts, sensitivity & IRFs)

The separate gamma diffuse dataset for training the energy regressor is only needed if you want to use the energy as feature in the particle classifier, but i'd recommend it.

You could add

• 2 gamma point-like datasets (optimize cuts, sensitivity & IRFs)

if you are interested only in point sources in the center of the field of view, where we have higher statistics thanks to the point like simulations.

[mdpunch]

So in the previous step, from six input merged files (3 gamma, 2 proton, 1 electron), I got one file each of {gamma/proton/electron}_final.dl2.h5, but for next step the procedure asks for two of each in input, for the optimization and IRF files respectively. So does that mean I have to use twelve files in input altogether to get two of each now? (since my assumption is that the files which are input here must have the ML applied).

Weren't you the one to remark that we shouldn't optimize the cuts on the same events as we compute the sensitivity? 😉

Indeed, but I was just "cheating" a bit here, to see if I could get the chain running at least!

Anyways, no, you only need one more of each, not twice the amount. You can also in principle, use the same files for each step, but as you correctly remarked, there is a danger of overfitting the cuts on the sensitivity.

So, if you want to use separate datasets for optimizing the cuts from computing sensitivity and IRFs, I think you will want:

* 4 gamma diffuse datasets (train energy regressor, train particle classifier, optimize cuts, sensitivity & IRFs)

* 3 proton datasets (train particle classifier, optimize cuts, sensitivity & IRFs)

* 2 electron datasets (optimize cuts, sensitivity & IRFs)

The separate gamma diffuse dataset for training the energy regressor is only needed if you want to use the energy as feature in the particle classifier, but i'd recommend it.

You could add

* 2 gamma point-like datasets (optimize cuts, sensitivity & IRFs)

if you are interested only in point sources in the center of the field of view, where we have higher statistics thanks to the point like simulations.

Thanks, I'll have to download a bunch more files then (planning on using ~10 files for each set, for these tests).

I am indeed only interested in point-like sources near the centre of the field, so does that mean I re-do the cut-optimization, sensitivity and IRFs with the point-like files, still using the energy regressor and particle classifier trained on the diffuse gammas?

[maxnoe]

Yes, training should always happen on diffuse, at least for the classifier, otherwise the classifier will mainly learn to distinguish point-like from diffuse, not gammas from protons.

[mdpunch]

So, if you want to use separate datasets for optimizing the cuts from computing sensitivity and IRFs, I think you will want:

* 4 gamma diffuse datasets (train energy regressor, train particle classifier, optimize cuts, sensitivity & IRFs)

* 3 proton datasets (train particle classifier, optimize cuts, sensitivity & IRFs)

* 2 electron datasets (optimize cuts, sensitivity & IRFs)

The separate gamma diffuse dataset for training the energy regressor is only needed if you want to use the energy as feature in the particle classifier, but i'd recommend it.
You could add

* 2 gamma point-like datasets (optimize cuts, sensitivity & IRFs)

if you are interested only in point sources in the center of the field of view, where we have higher statistics thanks to the point like simulations.

Coming back to this after a break (May holidays, SWAT tests, NectarCAM meeting organization, with a quick question.

Following the above, I have a bunch of merged files with obvious names:

gamma_diffuse_merged_train_en.std.dl2.h5
gamma_diffuse_merged_train_cls.std.dl2.h5
gamma_diffuse_merged_optimize_cuts.std.dl2.h5
gamma_diffuse_merged_irfs.std.dl2.h5
proton_merged_train_cls.std.dl2.h5
proton_merged_optimize_cuts.std.dl2.h5
proton_merged_irfs.std.dl2.h5
electron_merged_optimize_cuts.std.dl2.h5
electron_merged_irfs.std.dl2.h5
gamma_point_merged_optimize_cuts.std.dl2.h5
gamma_point_merged_irfs.std.dl2.h5

But for the first step, training the models as here: https://ctapipe.readthedocs.io/en/latest/user-guide/tools/dl2_guide.html there should be a test dataset, as well as the training dataset.

Is it legit to use the optimize_cuts datasets for the test part (since I guess the test outputs are just for showing the performance)?

[maxnoe]

But for the first step, training the models as here: https://ctapipe.readthedocs.io/en/latest/user-guide/tools/dl2_guide.html there should be a test dataset, as well as the training dataset.

Test here means performance evaluation. So in your case, I'd use the final data set for this, the one on which you compute IRFs, benchmarks, metrics etc.

[maxnoe]

This leads to my problem now, that there seems to be no ML information in the files (no RandomForest at least). See error message below.

Without knowing which commands you ran and with which configuration, it's hard to answer here.
After running ctapipe-apply-models with the trained models on your input files, the output files should have the columns corresponding to your configuration. The names can be set though and by default are taking from the model name.

Note that the output of the quickstart uses the scikit-learn ExtraTrees classifier and regressors, which are much faster to train than random forests. This is a configuration used for our tests to speed things up. For the actual analysis, we use random forests, you can find the config files here:
https://gitlab.cta-observatory.org/cta-computing/dpps/datapipe/pipeline-configurations/-/tree/main/ctapipe/v0.23/v1?ref_type=heads

[maxnoe]

Is there a train_disp_reconstructor.yaml corresponding to the RandomForest ?

No, but you can just replace all occurences of ExtraTreesRegressor_ with RandomForestRegressor_ in the config.

You can also change the estimator itself to a random forest for better performance and adjust depth / n_estimators like the other configs I linked.

[mdpunch]

Is there a train_disp_reconstructor.yaml corresponding to the RandomForest ?

No, but you can just replace all occurences of ExtraTreesRegressor_ with RandomForestRegressor_ in the config.

You can also change the estimator itself to a random forest for better performance and adjust depth / n_estimators like the other configs I linked.

It works to change to RandomForestRegressor_.

I think the estimator itself is also a RandomForest (I'm using your linked configs above).

[maxnoe]

For disp you used the quick start config, right?
In that config, norm_cls and sign_cls are the model types used for disp.

[maxnoe]

ctapipe/src/ctapipe/resources/train_disp_reconstructor.yaml

Lines 27 to 37 in 4e55c3e

	norm_cls: ExtraTreesRegressor
	norm_config:
	n_estimators: 10
	max_depth: 10
	log_target: True
	# All classification algorithms in scikit-learn are supported
	# (https://scikit-learn.org/stable/modules/classes.html)
	sign_cls: ExtraTreesClassifier
	sign_config:

[mdpunch]

Yes, that's what I changed to put RandomForest_ instead of ExtraTrees_

[mdpunch]

I have uploaded a repository to git on this:
https://github.com/mdpunch/ReducedWindow_IRFs/tree/main

For now, with just the few tens of files on my external disk, and using Max's @maxnoe performance script (with lots of help from him), the statistical fluctuations show that there aren't enough to get a proper sensitivity, but at least there is no big effect in the range where there are enough stats.
(compare the outputs in https://github.com/mdpunch/ReducedWindow_IRFs/blob/main/performance.ipynb and https://github.com/mdpunch/ReducedWindow_IRFs/blob/main/performance_ReducedWindow.ipynb)
e.g. sensitivities for standard window:

... and reduced windows:

Question for ASWG: Could we do a production of DL2 files with e.g. the pointing used for CTAO-N?