Error while trying to run lstchain_mc_r0_to_dl1 script to prod5 MC

[mari-taka]

I tried to run lstchain_mc_r0_to_dl1 script to prod5 MC that includes 4 LSTs and MAGIC, but I got the following error. I'm using lstchain 0.6.0 and ctapipe 0.8.0.

cp01 ~ $ lstchain_mc_r0_to_dl1 --input-file /fefs/aswg/workspace/MC_common/corsika7.7_simtelarray_20200629/prod5/4LSTs_MAGIC/gamma/zenith_20deg/south_pointing/sim_telarray/off0.4deg/data/gamma_20deg_180deg_run1___cta-prod5-lapalma_4LSTs_MAGIC_desert-2158m_4LSTs_MAGIC_mono_off0.4.simtel.gz

Traceback (most recent call last):
  File "/home/mari.takahasi/Work/miniconda3/envs/lst-dev/bin/lstchain_mc_r0_to_dl1", line 8, in <module>
    sys.exit(main())
  File "/home/mari.takahasi/Work/miniconda3/envs/lst-dev/lib/python3.7/site-packages/lstchain/scripts/lstchain_mc_r0_to_dl1.py", line 69, in main
    custom_config=config,
  File "/home/mari.takahasi/Work/miniconda3/envs/lst-dev/lib/python3.7/site-packages/lstchain/reco/r0_to_dl1.py", line 301, in r0_to_dl1
    write_array_info_08(subarray, output_filename)
  File "/home/mari.takahasi/Work/miniconda3/envs/lst-dev/lib/python3.7/site-packages/astropy/utils/decorators.py", line 124, in deprecated_func
    return func(*args, **kwargs)
  File "/home/mari.takahasi/Work/miniconda3/envs/lst-dev/lib/python3.7/site-packages/lstchain/io/io.py", line 368, in write_array_info_08
    stage1._write_instrument_configuration(subarray)
  File "/home/mari.takahasi/Work/miniconda3/envs/lst-dev/lib/python3.7/site-packages/ctapipe/tools/stage1.py", line 407, in _write_instrument_configuration
    serialize_meta=serialize_meta,
  File "/home/mari.takahasi/Work/miniconda3/envs/lst-dev/lib/python3.7/site-packages/astropy/table/connect.py", line 113, in __call__
    registry.write(instance, *args, **kwargs)
  File "/home/mari.takahasi/Work/miniconda3/envs/lst-dev/lib/python3.7/site-packages/astropy/io/registry.py", line 566, in write
    writer(data, *args, **kwargs)
  File "/home/mari.takahasi/Work/miniconda3/envs/lst-dev/lib/python3.7/site-packages/astropy/io/misc/hdf5.py", line 306, in write_table_hdf5
    serialize_meta=serialize_meta)
  File "/home/mari.takahasi/Work/miniconda3/envs/lst-dev/lib/python3.7/site-packages/astropy/io/misc/hdf5.py", line 323, in write_table_hdf5
    raise OSError(f"Table {path} already exists")
OSError: Table /configuration/instrument/telescope/camera/geometry_MAGICCam already exists

It seems the error happened while writing tables about array information to the DL1 file.

write_array_info(subarray, output_filename)
write_array_info_08(subarray, output_filename)

When writing /instrument/telescope/camera/{camera} in write_array_info(), it is skipped if {camera} already exists. Meanwhile, when writing /configuration/instrument/telescope/camera/{camera} in write_array_info_08(), it is not. "subarray.telescope_types" has two MAGICs, so I think MAGICCam was tried to be written twice and the error happened.

[TelescopeDescription(type=LST, name=MAGIC, optics=MAGIC, camera=MAGICCam),
 TelescopeDescription(type=LST, name=MAGIC, optics=MAGIC, camera=MAGICCam),
 TelescopeDescription(type=LST, name=LST, optics=LST, camera=LSTCam)]

I was able to solve it modifying write_array_info_08() like write_array_info().

def write_array_info_08(subarray, output_filename):
    """
    Write the array info to a ctapipe v0.8 compatible DL1 HDF5 file

    Parameters
    ----------
    subarray: `ctapipe.instrument.subarray.SubarrayDescription`
    output_filename: str
    """

    serialize_meta = True

    subarray.to_table().write(
      output_filename,
      path="/configuration/instrument/subarray/layout",
      serialize_meta=serialize_meta,
      append=True,
    )

    subarray.to_table(kind="optics").write(
      output_filename,
      path="/configuration/instrument/telescope/optics",
      append=True,
      serialize_meta=serialize_meta,
    )

    for telescope_type in subarray.telescope_types:
      ids = set(subarray.get_tel_ids_for_type(telescope_type))
      if len(ids) > 0: # only write if there is a telescope with this camera
        tel_id = list(ids)[0]
        camera = subarray.tel[tel_id].camera
        camera_name = f'geometry_{camera}'
        with tables.open_file(output_filename, mode='a') as f:
          telescope_chidren = f.root['/configuration/instrument/telescope']._v_children.keys()
          if 'camera' in telescope_chidren:
            cameras_name = f.root['/configuration/instrument/telescope/camera']._v_children.keys()
            if camera_name in cameras_name:
              print(
                f'WARNING during lstchain.io.write_array_info_08():',
                f'camera {camera_name} seems to be already present in the h5 file.'
              )
              continue
        camera.geometry.to_table().write(
          output_filename,
          path=f"/configuration/instrument/telescope/camera/geometry_{camera}",
          append=True,
          serialize_meta=serialize_meta          
        )
        camera.readout.to_table().write(
          output_filename,
          path=f"/configuration/instrument/telescope/camera/readout_{camera}",
          append=True,
          serialize_meta=serialize_meta
        )

What do you think? Best regards, Mari

[rlopezcoto]

Thank you very much @mari-taka for reporting a problem already with the proposed solution! As you are saying, this seems to be related to the MAGICCam addition to the new MC, not specifically any of the Prod5 parameters. Could you please open a PR with the proposed solution implemented in a similar way as it was done in write_array_info so it can be reviewed? Thanks

[maxnoe]

Can you tell for which file this happened? I suspect this is a bug also in ctapipe since the code was copied from there.

I think the error is actually assuming that two different telescope descriptions never have the same camera.

The two MAGICs here must be different, otherwise there would only be one entry in telescope_types.

[YoshikiOhtani]

Hi @MaxNoe, she tried to process the following data that I produced in the La Palma IT Container:

/fefs/aswg/workspace/MC_common/corsika7.7_simtelarray_20200629/prod5/4LSTs_MAGIC/gamma/zenith_20deg/south_pointing/sim_telarray/off0.4deg/data/gamma_20deg_180deg_run1___cta-prod5-lapalma_4LSTs_MAGIC_desert-2158m_4LSTs_MAGIC_mono_off0.4.simtel.gz

[maxnoe]

Why does a file called "4LSTs_MAGIC_mono" have two MAGIC telescopes?

[maxnoe]

Okay, the two MAGIC telescopes in that file are different from each other but have the same camera.

In [11]: s.subarray.tel[5] == s.subarray.tel[6]                                                                                               
Out[11]: False

In [12]: s.subarray.tel[5].camera == s.subarray.tel[6].camera                                                                                 
Out[12]: True

This is something not anticipated by the subarray writing code.

There is also another problem with this: two different telescopes should not share the same string representation, as this is used for the lookup.

We could and maybe should go to a tel_id based lookup.

I'll open issues over at ctapipe.

[maxnoe]

The difference is the mirror area:

In [14]: s.subarray.tel[5].optics                                                                                                             
Out[14]: OpticsDescription(name=MAGIC, equivalent_focal_length=16.97 m, num_mirros=1, mirror_area=231.94 m2)

In [15]: s.subarray.tel[6].optics                                                                                                             
Out[15]: OpticsDescription(name=MAGIC, equivalent_focal_length=16.97 m, num_mirros=1, mirror_area=235.22 m2)

[YoshikiOhtani]

Why does a file called "4LSTs_MAGIC_mono" have two MAGIC telescopes?

If I understood your question correctly, the answer is that the word "mono" means we allow the mono trigger to all telescopes, so it doesn't mean there is one MAGIC telescope. The reason we implemented the MAGIC telescope in the simulation is that, as far as I understand, we can (i) study the LST1-MAGIC observation performance, (ii) use the MCs to analyze the real data taken together with MAGIC.
Please correct me @rlopezcoto if I said something wrong...

[maxnoe]

Ah ok, I just assumed "MAGIC mono" means one MAGIC telescope, sorry, makes perfect sense.

[YoshikiOhtani]

And thank you very much for checking the data so quickly! So in my understanding, @mari-taka should not open the pull request and wait for the fixes on the ctapipe, because the two MAGIC description should be handled as the different telescopes, is it right?

[maxnoe]

No, because lstchain has copied the code here, the issue also directly affects lstchain.

I will make a fix in ctapipe, which can then be also copied here.

[YoshikiOhtani]

Ok, thanks!

[maxnoe]

Fixed here: https://github.com/cta-observatory/ctapipe/pull/1423

[maxnoe]

There is another issue related to this: It is not forseen by ctapipe, that two telescope descriptions have the same name at the same time.

But the two MAGIC descriptions look exactly the same, only the mirror area is slightly different.

Reading back the DL1 data will result in the wrong mirror area for one of the telescopes, since the lookup is made on the name of the telescope LST_MAGIC_MAGICCam, which is identical for both.

[rlopezcoto]

@YoshikiOhtani, @mari-taka, @MaxNoe anybody willing to copy the ctapipe solution into lstchain until a new ctapipe version is release?

[garciagenrique]

Fixed here: cta-observatory/ctapipe#1423

However, this will be implemented on ctapipe v0.9.1, right ? (could you confirm please ? @maxnoe).

So by implementing @mari-taka 's method (modifying the whole write_array_info_08 function) I was able to run the r0_to_dl1 stage for Prod5 Mcs. However, it also run (at gave the exactly same output), just by taking out the whole write_array_info_08 from r0_to_dl1.py.

If anybody could double check it, it will be nice. Then we could open a PR.

[maxnoe]

Yes, the PR with the fix is in ctapipe 0.9.1, released yesterday.

However, this does not fix this issue completely, since the matching of telescope description components when reading is still done by name, which here is not unique.

So when reading back the written subarray, one of the two different magics will be lost. As far as I know, the only thing that uses mirror area is the muon analysis, which happens before, so this might not be a real problem.

However, this issue will be fixed in the next ctapipe version by switching to an index based lookup and storage instead of this name based storage

[moralejo]

Hi @garciagenrique , if we go provisionally for @mari-taka 's solution we could process the needed MC, right? And both cameras would be there (although I agree with @maxnoe , probably mirror area is not used anywhere afterwards)

[garciagenrique]

Yes, we will be able to process them !

[rlopezcoto]

I also support exporting that part of the code into lstchain, otherwise moving to ctapipe v0.9.1 would be too much trouble. Whenever this PR is merged, we can make (yet another) bugfix release to continue with the MC processing.

There is also a small issue with astropy incompatibilities with those from astropy. @garciagenrique and I worked it out yesterday night and he will also open a PR with the solution

[garciagenrique]

@mari-taka, do you want to create the PR to fix this ? Otherwise I/we can do it.

[maxnoe]

Actually, the change to ctapipe 0.9 should be rather seemless, nothing compared to the 0.7 -> 0.8 transition

[mari-taka]

Hi, @garciagenrique, I can create the PR, if it is fine for you. I think it does not take much time as I already have a solution.

[maxnoe]

We now rely on Subarray.{to,from}_hdf to read and write subarray information. The issue with the MAGIC telescopes will be fixed in the new data model when switching to ctapipe 0.12