bump version of ctapipe-extra directory to v0.3.2

[kosack]

This is the same as v0.3.1, but with camera readout definitions from prod3b

(the data have been updated on the dataserver)

[kosack]

Partially solves #1450, and will also help #1568

[maxnoe]

@kosack something seems to have changed for the LSTCam that makes tests fail.

[kosack]

Strangely the test for concentration passes on my machine. But I do get a failure for test_data_volume_reducer.

[kosack]

The only changes in the new files are:

• the camreadout FITS files exist (before they did not, so a dummy signal pulse was used)
• the camgeom fits files seem to be identical, except the
  • pixel area has changed (not sure why - do we have a bug?) it changed from 0.00216 deg2 to 0.00207 deg2. That change will affect the toymodel, since it scales the intensity by pixel area (so the intensities will reduce a bit)
  • the PIX_TYPE header changed from "hexagonal" to "hexagon" (need to check if this causes any problems anywhere) [edit: I don't see how it could, it's now an ENUM]

The reference pulse is used during image extraction, if the integration correction is enabled (normally yes). So for tests that did not load up a real event, perhaps this causes some difference?
I do note that some tests are only ever tested against the LSTCam, which could mean they are hard-coded to work only there and should be expanded to test other camers (e.g. test_data_volume_reducer).

[kosack]

The problem with test_data_volume_reducer is that the test assumes that the signal will be a strict sum of the input waveform, but since now with the reference pulse, the integration correction is performed, the value is not exactly that anymore, so it fails. Even with apply_integration_correction turned off, however, I get a different answer than expected.

The following config seems to fix it (the tests assumes a really simple integration, so I use FixedWindowSum):

 {
            "TailCutsDataVolumeReducer": {
                "TailcutsImageCleaner": {
                    "picture_threshold_pe": 700.0,
                    "boundary_threshold_pe": 350.0,
                    "min_picture_neighbors": 0,
                    "keep_isolated_pixels": True,
                },
                "image_extractor_type": "FixedWindowSum",
                "FixedWindowSum": {
                    "apply_integration_correction": False,
                },
                "n_end_dilates": 1,
                "do_boundary_dilation": True,
            }
        }

@Hckjs is this ok?

[jsitarek]

Hi @kosack
I was looking into this in parallel and also suspected the integration correction.
In the fix above you changed the
NeighborPeakWindowSum to FixedWindowSum extractor
the reason why it was not working with the first extractor is that the waveforms in the test are flat, [100, 100, ... 100] so it selects as the peak position 0th sample. The integration is with 7 slices (hence the code expects 700), however the default of window_shift is 3, so it would integrate from -3 up to 3rd slice, and the negative ones are skipped so you end up with a sum of 400.

alternative solution if "NeighborPeakWindowSum" is to be used is to set:
"NeighborPeakWindowSum": {
"apply_integration_correction": False,
"window_shift":0

this worked for me

[kosack]

Thanks, yes, we both notice the same thing - in fact I like your fix better, since it still uses the default extractor

[jsitarek]

should I commit it (to the other branch where we had the same problem?)
this still leaves the problem with the concentration test

[jsitarek]

I also have a suspicion about the problem with concentration. When I run the test on my PC it gives the value of 0.052 with the limits 0.05-0.2, while in the automatic test it gave 0.0477 just below the limit.
The test is using create_sample_image, which is using generate_image from toymodel.py and there there are
signal = np.random.poisson(expected_signal)
noise = np.random.poisson(nsb_level_pe, size=signal.shape)

I'm not sure if the default random generator is set to any specific seed for all the test, but if no this could explain the problem any why it is not reproducable. The cleanest solution I guess would be to define a random number generator inside the toymodel class and pass the seed to it before the test. However there might be more tests that depend on the same image...

[maxnoe]

The cleanest solution I guess would be to define a random number generator inside the toymodel class and pass the seed to it before the test

Just add np.random.seed(<your favourite seed>) at the top of the test function.

[jsitarek]

sorry, I was actually mistaken, the seed is set inside the create_sample_image, so it must be something different, unless something funny happens with the random number generator (e.g. that the same one is used in parallel to different tests)

[maxnoe]

@jsitarek The tests do not run in parallel.

I think the change in pixel area is the culprit. That explains the slightly lower value.

[jsitarek]

Hi @maxnoe
if this is the area, why a different values os obtained when testing locally and when running the automatic checks?

[maxnoe]

@jsitarek I guess because the old file is still in the cache?

rm -rf $HOME/.cache/ctapipe

The cache only works by filename.

[kosack]

I'll commit it here - it's best to keep it separate.

[kosack]

The cache only works by filename.

For safety, maybe we should check more than filename, like a checksum.

[kosack]

I think the change in pixel area is the culprit. That explains the slightly lower value.

Yes, this seems to be the issue - the question is why it changed, when the pixel positions are the same. The algorithm implementation to compute it was refactored, but I checked the version used to make the original, and the formulas are identical, just moved around a bit.

[jsitarek]

sorry @maxnoe but I do not understand the possible problem with the pixel sizes
I cleaned the cache in my PC, however when I rerun the test I still get the same value 0.052 that passes the test. When the test is run automatically I guess cache is automatically cleaned

EDIT: also the concentration that was triggering the difference is just defined as the brightest pixel to the total light intensity, and this does not depend on the area. What am I missing?

[maxnoe]

For safety, maybe we should check more than filename, like a checksum.

That would require always downloading at least the catalog of checksums. Currently the download is only implemented as fallback.

Simpler would be to store the cache by url and check that the url is the same.

[kosack]

Pix area is the main difference:

curl http://cccta-dataserver.in2p3.fr/data/ctapipe-extra/v0.3.1/LSTCam.camgeom.fits.gz -o LSTCam031.camgeom.fits.gz
curl http://cccta-dataserver.in2p3.fr/data/ctapipe-extra/v0.3.2/LSTCam.camgeom.fits.gz -o LSTCam032.camgeom.fits.gz
fitsdiff LSTCam031.camgeom.fits.gz LSTCam032.camgeom.fits.gz

fitsdiff LSTCam031.camgeom.fits.gz LSTCam032.camgeom.fits.gz -a 0.00001

 fitsdiff: 4.2
 a: LSTCam031.camgeom.fits.gz
 b: LSTCam032.camgeom.fits.gz
 Maximum number of different data values to be reported: 10
 Relative tolerance: 0.0, Absolute tolerance: 1e-05

Extension HDU 1:

   Headers contain differences:
     Keyword PIX_TYPE has different values:
        a> hexagonal
         ?        --
        b> hexagon
     Keyword SOURCE   has different values:
        a> gamma_20deg_180deg_run1000___cta-prod3-demo_desert-2150m-Paranal-demo2rad_cone10.simtel.gz
        b> /Users/kkosack/Data/CTA/Prod3b/proton/proton_20deg_0deg_run10___cta-prod3_desert-2150m-Paranal-merged.simtel.gz
     Keyword TAB_VER  has different values:
        a> 1.0
        b> 1.1

   Data contains differences:
     Column pix_area data differs in row 0:
        a> 0.0021649681944033382
        b> 0.002079326892271638
     Column pix_area data differs in row 1:
        a> 0.0021649681944033382
        b> 0.002079326892271638
     Column pix_area data differs in row 2:
        a> 0.0021649681944033382
        b> 0.002079326892271638
     Column pix_area data differs in row 3:
        a> 0.0021649681944033382
        b> 0.002079326892271638
     Column pix_area data differs in row 4:
        a> 0.0021649681944033382
        b> 0.002079326892271638
     Column pix_area data differs in row 5:
        a> 0.0021649681944033382
        b> 0.002079326892271638
     Column pix_area data differs in row 6:
        a> 0.0021649681944033382
        b> 0.002079326892271638
     Column pix_area data differs in row 7:
        a> 0.0021649681944033382
        b> 0.002079326892271638
     Column pix_area data differs in row 8:
        a> 0.0021649681944033382
        b> 0.002079326892271638
     Column pix_area data differs in row 9:
        a> 0.0021649681944033382
        b> 0.002079326892271638
     ...1845 additional difference(s) found.
     ...
     1855 different table data element(s) found (25.00% different).

[kosack]

So if we think the new pixel areas are correct, we should just update the test values of the concentration test - I don't think the randomness matters too much, it should be within tolerance. But I would like to understand why the area is not the same. A 4% change is quite strange.

[maxnoe]

@kosack Could it be that that simply changed between prod2 and prod3?

The older data comes from gamma_test.simtel.gz and In think that was prod2.

[kosack]

The area is not stored in the files, it's computed from the pixel positions (and those are the same between the two files)

[maxnoe]

@kosack that is not true, the area is read from the files.

ctapipe/ctapipe/io/simteleventsource.py

Lines 65 to 88 in b5c90a1

	def build_camera(cam_settings, pixel_settings, telescope, frame):
	pixel_shape = cam_settings["pixel_shape"][0]
	try:
	pix_type, pix_rotation = CameraGeometry.simtel_shape_to_type(pixel_shape)
	except ValueError:
	warnings.warn(
	f"Unkown pixel_shape {pixel_shape} for camera_type {telescope.camera_name}",
	UnknownPixelShapeWarning,
	)
	pix_type = "hexagon"
	pix_rotation = "0d"
	geometry = CameraGeometry(
	telescope.camera_name,
	pix_id=np.arange(cam_settings["n_pixels"]),
	pix_x=u.Quantity(cam_settings["pixel_x"], u.m),
	pix_y=u.Quantity(cam_settings["pixel_y"], u.m),
	pix_area=u.Quantity(cam_settings["pixel_area"], u.m ** 2),
	pix_type=pix_type,
	pix_rotation=pix_rotation,
	cam_rotation=-Angle(cam_settings["cam_rot"], u.rad),
	apply_derotation=True,
	frame=frame,
	)

[kosack]

@kosack that is not true, the area is read from the files.

Aha, then that must be the difference: before it was guessed from the pixel distances, but of course in reality the area is not exactly that since there is some width to the winston cone. That answers the question: I haven' realized we now read it (though it makes sense now that we use pyeventio). So the new areas are more correct. Therefore, let's just update the test values. In fact, it's not a very good test. It would be better to make an artificial camera image (not a random toy model)

[kosack]

I just changed the test values. I think we should open an issue to create a better test for that, using an artificial input where the concentration is known a priori, but that can come later.

[kosack]

Now there is a problem with toymodel.test_intensity, which seems quite strange:

 assert np.average(geom.pix_x.to_value(u.m), weights=signal) == approx(0.2, rel=0.15)

 0.23126951261659126 == 0.2 ± 3.0e-02

This test checks that a gaussian image centered at 0.2 along the x-axis has the intensity values one would expect. The only thing that has changed here is the pixel area (3% lower than before), but that should only affect the overall scale of the image, not the first-order moments of the distribution.

However, I think the issue is a very low image intensity. The "signal" value here is :

>>> print(signal[signal>0])
array([1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 4, 2, 1, 2, 1, 1,1, 1, 1, 2, 1, 1, 1, 1, 1]) 

So clearly not very intense! Increasing the total intensity of the model to 200 PE, it works fine. So I'll do that, I think the test was only passing by chance before. @maxnoe any objection?

[maxnoe]

@kosack no, looks good

[codecov]

Codecov Report

Merging #1588 (9e6ae3e) into master (b5c90a1) will decrease coverage by 0.02%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master    #1588      +/-   ##
==========================================
- Coverage   90.81%   90.78%   -0.03%     
==========================================
  Files         192      192              
  Lines       13966    13966              
==========================================
- Hits        12683    12679       -4     
- Misses       1283     1287       +4     

Impacted Files	Coverage Δ
ctapipe/image/tests/test_reducer.py	100.00% <ø> (ø)
ctapipe/image/tests/test_concentration.py	95.23% <100.00%> (ø)
ctapipe/image/tests/test_toy.py	100.00% <100.00%> (ø)
ctapipe/utils/datasets.py	89.10% <100.00%> (-1.00%)	⬇️
ctapipe/instrument/camera/readout.py	85.50% <0.00%> (-4.35%)	⬇️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update b5c90a1...ce1be58. Read the comment docs.

[maxnoe]

@kosack Just to confirm:

In [9]: cam.pix_area
Out[9]: 
<Quantity [0.00207933, 0.00207933, 0.00207933, ..., 0.00207933, 0.00207933,
           0.00207933] m2>

In [10]: cam.guess_pixel_area(cam.pix_x, cam.pix_y, cam.pix_type)
Out[10]: 
<Quantity [0.00216497, 0.00216497, 0.00216497, ..., 0.00216497, 0.00216497,
           0.00216497] m2>

[kosack]

Yes, the guess method doesn't perfectly reproduce the real areas, since in reality some intensity is lost to gaps/edges of the winston cone. We could add a fudge factor for that to make it more realistic, by comparing the real values to the guessed values for all cameras, and averaging for example. In that case we'd just multiply the guessed area by e.g. 98%. But that should be a separate issue/PR.

The areas are mostly only used for plotting, so the effect can also be ignored. The only place so far they are used for computation is in the toy model, but that doesn't need to be perfect.

[maxnoe]

The areas are mostly only used for plotting, so the effect can also be ignored. The only place so far they are used for computation is in the toy model, but that doesn't need to be perfect.

Muon intensity fitter.