Error running lstchain_mc_sensitivity.py script

[andres-baquero]

I'm using lstchain v0.4.4 and I've installed ctaplot.

I'm getting this error running lstchain_mc_sensitivity.py:

`AttributeErrorTraceback (most recent call last)
in
6 n_bins_energy, n_bins_gammaness,
7 n_bins_theta2, noff,
----> 8 obstime)

/fefs/home/andres.baquero/lstchain-dev/lstchain_v044/cta-lstchain/lstchain/mc/sensitivity.py in find_best_cuts_sensitivity(simtelfile_gammas, simtelfile_protons, dl2_file_g, dl2_file_p, nfiles_gammas, nfiles_protons, n_bins_energy, n_bins_gammaness, n_bins_theta2, noff, obstime)
332
333 # Pass units to GeV and cm2
--> 334 mc_par_g['emin'] = mc_par_g['emin'].to(u.GeV)
335 mc_par_g['emax'] = mc_par_g['emax'].to(u.GeV)
336

AttributeError: 'numpy.float64' object has no attribute 'to'
`
Thanks for your help

[rlopezcoto]

Hi @andres-baquero, the minimum energy (mc_par_g['emin']) and the maximum energy (mc_par_g['emax']) should be read from the original MC file (simtelfile_gammas), which one is the one you are using?

[andres-baquero]

Hi @rlopezcoto. I'm using a dl1 file of gamma point merged with all runs produced from:
/fefs/aswg/data/mc/DL0/20190415/gamma/south_pointing

[rlopezcoto]

Those are DL0 according to the directory, right?

[andres-baquero]

Yes. I'm using DL1 files for protons and gammas that have been produced from this DL0 files.

[rlopezcoto]

Can you give me the exact command you are using (with all the input files) so I can try to reproduce the error?
Thanks

[andres-baquero]

Thanks @rlopezcoto for your help. I'm using this:

python /fefs/aswg/software/virtual_env/ctasoft/cta-lstchain/lstchain/scripts/lstchain_mc_sensitivity.py -gd1 /fefs/aswg/data-test/mc/20190415/dl1/v0.4.3_v00/gamma/tailcut_6_3/South_pointing/merge/dl1_gamma_6_3_merged_mc_tables.h5 -pd1 /fefs/aswg/data-test/mc/20190415/dl1/v0.4.3_v00/proton/tailcut_6_3/South_pointing/merge/dl1_proton_6_3_merged_mc_tables.h5 -gd2-cuts /fefs/aswg/data-test/mc/20190415/dl2/irf_test/dl2_dl1_gamma_6_3_merged_mc_tables.h5 -pd2-cuts /fefs/aswg/data-test/mc/20190415/dl2/irf_test/dl2_dl1_proton_6_3_merged_mc_tables.h5 -gd2-sens /fefs/aswg/data-test/mc/20190415/dl2/irf_test/dl2_dl1_gamma_6_3_merged_mc_tables_sens.h5 -pd2-sens /fefs/aswg/data-test/mc/20190415/dl2/irf_test/dl2_dl1_proton_6_3_merged_mc_tables_sens.h5

[moralejo]

/fefs/home/andres.baquero/lstchain-dev/lstchain_v044/cta-lstchain/lstchain/mc/sensitivity.py in find_best_cuts_sensitivity(simtelfile_gammas, simtelfile_protons, dl2_file_g, dl2_file_p, nfiles_gammas, nfiles_protons, n_bins_energy, n_bins_gammaness, n_bins_theta2, noff, obstime)
332
333 # Pass units to GeV and cm2
--> 334 mc_par_g['emin'] = mc_par_g['emin'].to(u.GeV)
335 mc_par_g['emax'] = mc_par_g['emax'].to(u.GeV)
336

AttributeError: 'numpy.float64' object has no attribute 'to'

Did that ever work? Seems to me the values mc_par_g['emin'] are those directly read from the hdf5 file, so they have no associated units. Hence we should do mc_par_g['emin']*u.TeV (if the number stored in the file is in TeV), and then do the conversion to GeV.

[andres-baquero]

Hi @moralejo yes that was the problem in this case.

[vuillaut]

Hi @andres-baquero
If you fixed the bug, can you make a PR, please?
ps: please get everything in TeV as the default unit.

Thanks.

[andres-baquero]

Hi @vuillaut . Ok I will make the PR.

Just one question @rlopezcoto the script receive 4 DL2 files:

--gammadl2-cuts DL2_FILE_G_CUTS, --protondl2-cuts DL2_FILE_P_CUTS, --gammadl2-sens DL2_FILE_G_SENS, --protondl2-sens DL2_FILE_P_SENS

Which is the difference between this DL2 files ?.

[SeiyaNozaki]

Hi @andres-baquero ,
I reproduced the issue like below;

>>> from lstchain.mc.sensitivity import process_mc
>>> dl1_gamma="/fefs/aswg/data/mc/DL1/20190415/gamma/south_pointing/20191128_v.0.3.1_v00/dl1_gamma_south_pointing_20191128_v.0.3.1_v00_DL1_testing.h5"
>>> dl2_gamma="/fefs/aswg/data/mc/DL2/20190415/gamma/south_pointing/20191128_v.0.3.1_v00/dl2_dl1_gamma_south_pointing_20191128_v.0.3.1_v00_DL1_testing.h5"
>>> gammaness_g, theta2_g, e_reco_g, e_true_g, mc_par_g, events_g = process_mc(dl1_gamma, dl2_gamma, 'gamma')
Table '/simulation/run_config' has column 'obs_id' that is not in container MCHeaderContainer. It will be skipped
Table '/simulation/run_config' is missing column 'prod_site_array' that is in container MCHeaderContainer. It will be skipped.
Table '/simulation/run_config' is missing column 'prod_site_coord' that is in container MCHeaderContainer. It will be skipped.
Table '/simulation/run_config' is missing column 'prod_site_subarray' that is in container MCHeaderContainer. It will be skipped.
>>> mc_par_g['emin']
0.004999999888241291

But with DL1 file which I generated by using the latest lstchain, mc_par_g['emin'] has TeV unit.
I'm checking...

>>> dl1_gamma="/fefs/aswg/workspace/seiya.nozaki/data/MC/IC6_3/DL1/gamma_point/dl1_gamma_20deg_0deg_run1___cta-prod3-demo-2147m-LaPalma-baseline-mono_off0.4.simtel.h5"
>>> gammaness_g, theta2_g, e_reco_g, e_true_g, mc_par_g, events_g = process_mc(dl1_gamma, dl2_gamma, 'gamma')
Table '/simulation/run_config' has column 'obs_id' that is not in container MCHeaderContainer. It will be skipped
Table '/simulation/run_config' is missing column 'prod_site_array' that is in container MCHeaderContainer. It will be skipped.
Table '/simulation/run_config' is missing column 'prod_site_coord' that is in container MCHeaderContainer. It will be skipped.
Table '/simulation/run_config' is missing column 'prod_site_subarray' that is in container MCHeaderContainer. It will be skipped.
>>> mc_par_g['emin']
<Quantity 0.005 TeV>

[vuillaut]

Actually I am puzzled with the bug here.
Checking the code, mc_par_g['emin'] must have unit here.
It uses read_simu_info_hdf5 that returns a MCHeaderContainer() with TeV unit for energy_range_min. And this part of the code is acutally unit tested test_read_sim_par.

As pointed by @SeiyaNozaki, the issue here could be that the files were generated with an older version of lstchain where units were not applied (but in this case I would have expected the reading to fail - so before the attempted conversion).

As @moralejo also asked for it, we'll reduce the MC production with the latest release.

[andres-baquero]

So @vuillaut the PR is still needed ?.

[vuillaut]

So @vuillaut the PR is still needed ?.

If I am correct, actually not.

[vuillaut]

Hi @andres-baquero
I finally found the issue. It actually comes from the merging, which is done in a generic way for all tables, loosing the unit information. Nasty one :-/

[rlopezcoto]

Sorry for arriving so late at the discussion.
@moralejo: the hdf5 can have units, but as reported by @vuillaut, they were stripped away in the merging.
@andres-baquero: there are not 4 DL2 files, there are two sets of DL2 files, the ones to calculate the optimal theta2 and gammaness cuts (gammadl2-cuts protondl2-cuts) and the ones to apply them to calculate the sensitivity (gammadl2-sens protondl2-sens). This was set like this by @misabelber not to have a biased sensitivity.
@SeiyaNozaki: thanks for the checks!
@vuillaut: thanks a lot for spotting the problem, is #300 working it out or do we need another fix?

[rlopezcoto]

Fixed in #299