Problems with some subruns from 2023-08-19 in muon analysis

[morcuended]

There were some subruns from 2023-08-19 runs: 14034.0021 (around the middle of the run) and other 6 subruns of run 14029 (25, 26, 27, 29, 30, 31 towards the end of the run) that failed the R0 to DL1 step. Some problems (related to the rates) were reported in the ELOG during those runs.

Example of R0 to DL1 command to reproduce the error:

lstchain_data_r0_to_dl1 --input-file=/fefs/aswg/data/real/R0/20230819/LST-1.1.Run14034.0021.fits.fz --output-dir=/fefs/aswg/data/real/running_analysis/20230819/v0.10 --pedestal-file=/fefs/aswg/data/real/monitoring/PixelCalibration/Cat-A/drs4_baseline/20230819/pr
o/drs4_pedestal.Run14024.0000.h5 --calibration-file=/fefs/aswg/data/real/monitoring/PixelCalibration/Cat-A/calibration/20230819/v0.10.2/calibration_filters_52.Run14031.0000.h5 --time-calibration-file=/fefs/aswg/data/real/monitoring/PixelCalibration/Cat-A/drs4_time_samplin
g_from_FF/20220518/v0.10.0/time_calibration.Run08349.0000.h5 --systematic-correction-file=/fefs/aswg/data/real/monitoring/PixelCalibration/Cat-A/ffactor_systematics/20230410/v0.10.3/calibration_scan_fit_20230410.0000.h5 --config=/fefs/aswg/lstosa/config_files_lstchain/lst
chain_standard_v0.10_heuristic_ff.json --pointing-file=/fefs/onsite/monitoring/driveLST1/DrivePositioning/DrivePosition_log_20230819.txt --run-summary-path=/fefs/aswg/data/real/monitoring/RunSummary/RunSummary_20230819.ecsv

Here is the error traceback:

...
Traceback (most recent call last):
  File "/fefs/aswg/software/conda/envs/osa-v0.10/lib/python3.11/site-packages/lstchain/scripts/lstchain_data_r0_to_dl1.py", line 227, in main
    r0_to_dl1.r0_to_dl1(
  File "/fefs/aswg/software/conda/envs/osa-v0.10/lib/python3.11/site-packages/lstchain/reco/r0_to_dl1.py", line 701, in r0_to_dl1
    analyze_muon_event(subarray,
  File "/fefs/aswg/software/conda/envs/osa-v0.10/lib/python3.11/site-packages/lstchain/image/muon/muon_analysis.py", line 210, in analyze_muon_event
    muonringparam, clean_mask, dist, image_clean = fit_muon(x, y, image, geom,
                                                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/fefs/aswg/software/conda/envs/osa-v0.10/lib/python3.11/site-packages/lstchain/image/muon/muon_analysis.py", line 134, in fit_muon
    ring = fitter(x, y, image, clean_mask)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/fefs/aswg/software/conda/envs/osa-v0.10/lib/python3.11/site-packages/ctapipe/image/muon/ring_fitter.py", line 43, in __call__
    radius, center_fov_lon, center_fov_lat = fit_function(
                                             ^^^^^^^^^^^^^
  File "/fefs/aswg/software/conda/envs/osa-v0.10/lib/python3.11/site-packages/ctapipe/image/muon/ring_fitter.py", line 17, in kundu_chaudhuri
    return kundu_chaudhuri_circle_fit(fov_lon[mask], fov_lat[mask], weights[mask])
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/fefs/aswg/software/conda/envs/osa-v0.10/lib/python3.11/site-packages/ctapipe/image/muon/fitting.py", line 45, in kundu_chaudhuri_circle_fit
    np.sum(weights * ((center_x - x) ** 2 + (center_y - y) ** 2)) / weights_sum
                       ~~~~~~~~~^~~
  File "/fefs/aswg/software/conda/envs/osa-v0.10/lib/python3.11/site-packages/astropy/coordinates/angles.py", line 695, in __array_ufunc__
    results = super().__array_ufunc__(*args, **kwargs)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/fefs/aswg/software/conda/envs/osa-v0.10/lib/python3.11/site-packages/astropy/units/quantity.py", line 691, in __array_ufunc__
    raise e
  File "/fefs/aswg/software/conda/envs/osa-v0.10/lib/python3.11/site-packages/astropy/units/quantity.py", line 636, in __array_ufunc__
    converters, unit = converters_and_unit(function, method, *inputs)
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/fefs/aswg/software/conda/envs/osa-v0.10/lib/python3.11/site-packages/astropy/units/quantity_helper/converters.py", line 199, in converters_and_unit
    raise UnitConversionError(
astropy.units.core.UnitConversionError: Can only apply 'subtract' function to dimensionless quantities when other argument is not a quantity (unless the latter is all zero/infinity/nan).

It is the first time I see this error in the onsite data processing.

It made the script fail. Maybe it would be good to not fully exit the script due to the muon analysis and keep going with the rest of the steps. Although, it could indicate actual problems in the data taking.

I leave this issue here for the record.

[maxnoe]

Maybe it would be good to not fully exit the script due to the muon analysis and keep going with the rest of the steps. Although, it could indicate actual problems in the data taking.

Since this is the first time this happens and it seems to depend on specific data runs, I'd rather have this raise an exception and we investigate immediately rather than an error being catched, and we only see a slight drop in muon rate.

[maxnoe]

That being said, I am quite puzzled what problem in the data could result in a wrong or missing unit only in the muon code.

[maxnoe]

The issue is that both center_x and center_y don't have a unit, while x and y do not, it seems only for this event.

Investigating further, it seems that only for this event, weights is a masked array and the unit gets lost if that is the case:

ipdb> np.sum(weights.data * (x - mean_x) * x)
<Quantity 1171.52607276 deg2>
ipdb> np.sum(weights * (x - mean_x) * x)
1171.5260727590803

So the question is: why is weights suddenly a masked array?

[maxnoe]

weights is the dl1 image, image and clean_mask are both masked arrays for this event

[maxnoe]

The reason seems to be that event.mon.tel[telescope_id].calibration.unusable_pixels is a masked array

@FrancaCassol

[moralejo]

The reason seems to be that event.mon.tel[telescope_id].calibration.unusable_pixels is a masked array

@FrancaCassol

So is this related to the numpy issue you recently found? (and appears only now because of a change in numpy version)

[FrancaCassol]

no, I don't think, I am on it

[FrancaCassol]

The problem is due to these lines:

cta-lstchain/lstchain/reco/r0_to_dl1.py

Line 664 in 370435f

bad_pixels = event.mon.tel[telescope_id].calibration.unusable_pixels[0]

and

cta-lstchain/lstchain/reco/r0_to_dl1.py

Line 673 in 370435f

bad_pixels_hg = event.mon.tel[telescope_id].calibration.unusable_pixels[0]

cta-lstchain/lstchain/reco/r0_to_dl1.py

Line 674 in 370435f

bad_pixels_lg = event.mon.tel[telescope_id].calibration.unusable_pixels[1]

which use the event.mon.tel[telescope_id] to calibrate the muons, but this is updated by the interleave data processing when enough calibration statistics is reached. In general, this happens at the end of the sub_run because the sub_run interleave statistics is always less than the requested 10000 events, however for this run it seems that the 10000 statistics is reached (?). Then, the unusable_pixel array is updated and it becomes a masked array, giving the problem. This also must be corrected. I will open a PR with all the corrections.

[moralejo]

Is it possible that it reaches the number just because of wrong tagging, i.e. those 10000 are not really interleaved events?

[FrancaCassol]

it should not be, because those events are tagged as UNKNOWN, but we must verify it

[FrancaCassol]

There are indeed 10134 FF events in this subrun (trigger_type 1), and there are 823 events with trigger_type 4 (calibration trigger), tagged with event_type 255 (unknown events because they do not fit the heuristic criteria)

[morcuended]

fixed in #1151