Refactoring coreas to include interpolation #688

cg-laser · 2024-06-05T13:12:03Z

Major refactoring of the coreas interface to read the electric field from coreas hdf5 files.
A new class readCoREASDetector.py can be used to read an array and get the interpolated position.

This is largely a copy of #685 which is an undefined state because of an accidental force push. Please never force-push unless you know what you are doing.

…ing readCoREAS.py which reads a random position and readCoREASStationGrid.py which takes the closest observer in an array. Now readCoREASDetector.py can be used to read an array and get the interpolated position.

…na position) code seems to give expected results now.

… improve docstring

cg-laser · 2024-06-05T18:47:41Z

@MijnheerD @lpyras I'm getting these warning from the interpolator module all the time:

Trace arrival times were not set during init, only relative timings are returned!
warning: negative values in abs_spectrum found: 6 times. Setting to zero.

are we not using the interpolator correctly?

…test will be part of PR #689

cg-laser · 2024-06-14T13:32:45Z

@MijnheerD @lpyras I added an example of a "LOFAR"-style Xmax reco. The example is using a vertical shower. It would probably be good to also test with an inclined shower as some potential bugs might only appear then.

Before merging this PR, it would be good if the module would be tested a bit more. Is my understanding correct, that this module has not been used by anyone and is essentially untested?

cg-laser · 2024-06-14T13:34:32Z

The output of the reconstruction looks like this.
A remaining ToDo is a correct uncertainty estimation of the energy fluence. This should be implemented into the efieldSignalReconstructor

…ation.

lpyras · 2024-08-05T22:26:51Z

@MijnheerD @lpyras I'm getting these warning from the interpolator module all the time:
Trace arrival times were not set during init, only relative timings are returned!
warning: negative values in abs_spectrum found: 6 times. Setting to zero.
are we not using the interpolator correctly?

@MijnheerD I dont really know what that means. I always thought we are anyway only interested in the relative timing. I implemented access now the timing information provided by the cr-pulse-interpolator. Is there something else to do?

lpyras · 2024-08-05T22:27:44Z

@MijnheerD @lpyras I added an example of a "LOFAR"-style Xmax reco. The example is using a vertical shower. It would probably be good to also test with an inclined shower as some potential bugs might only appear then.

Before merging this PR, it would be good if the module would be tested a bit more. Is my understanding correct, that this module has not been used by anyone and is essentially untested?

I tested it before the coreas refactoring for effective area calculation. That yielded reasonable results.

lpyras · 2024-08-06T02:12:36Z

I changed everything we discussed so far e.g. change the coreas readers in a way that they only rely on the readCorsika7() function. Unfortunately I didn't have time to test everything. Any volunteers? Otherwise I try to do it in the next weeks. :)

MijnheerD · 2024-08-07T09:47:25Z

Sorry for coming in late on this, I thought I already answered some the questions but apparently I didn't.

@MijnheerD @lpyras I'm getting these warning from the interpolator module all the time:
Trace arrival times were not set during init, only relative timings are returned!
warning: negative values in abs_spectrum found: 6 times. Setting to zero.
are we not using the interpolator correctly?
@MijnheerD I dont really know what that means. I always thought we are anyway only interested in the relative timing. I implemented access now the timing information provided by the cr-pulse-interpolator. Is there something else to do?

This is important, because usually the relative timing is not enough. The problem is that for larger starshapes, the CoREAS time traces start at different times. We added the option to pass the start times to the interpolator, which it can then interpolate for you as well. The warning tells you that the start times were not provided to the constructor. On the lofar/interpolator branch we added this on line 65 of "readCoREASInterpolator.py". Then on line 191 we retrieve the start times, which are used in the make_sim_station function give the traces the correct start time.

@MijnheerD @lpyras I added an example of a "LOFAR"-style Xmax reco. The example is using a vertical shower. It would probably be good to also test with an inclined shower as some potential bugs might only appear then.

Before merging this PR, it would be good if the module would be tested a bit more. Is my understanding correct, that this module has not been used by anyone and is essentially untested?

I agree we should test it more, and I think once the LOFAR pipeline is done we will have an excellent playground for that. Although we can also already start testing on simulations, like the Xmax example (thanks @cg-laser for the script!). I am not sure when we will be able to get to it, but in September we should have a master student who will need to work with the interpolator so he might be able to test it more thoroughly.

anelles · 2024-08-07T14:44:33Z

I just tested the branch and it actually fails with a pip installed cr_interpolator on:
File "/usr/local/lib/python3.12/site-packages/cr_pulse_interpolator/signal_interpolation_fourier.py", line 9, in
import interpolation_fourier as interpF
ModuleNotFoundError: No module named 'interpolation_fourier
I feel like I have seen this error discussed. Do we need a new cr_interpolator release?

fschlueter · 2025-01-29T11:42:16Z

I am happy to approve this PR. I did not go line by line through this monolithic PR but I hope if people start using it, necessary fixes will be made. I have one question though. Besides adding the core functionality of the interpolator module this PR also removes examples. Why is that?

MijnheerD · 2025-01-29T12:24:50Z

The change I just made should hopefully fix the last unit test 🙏🏽

MijnheerD · 2025-01-29T12:26:47Z

With this PR we will get rid of the unit test "Tiny reconstrucution". Should we add a new test with the file? And should we do this in this PR?

Also, do we still need the example SimpleMCReconstruction.py and CostumHybridDetector.py? Both break with the current changes. I made new examples in cr_analysis. But maybe they are hard to find? Should I call them air_shower_scripts?

The SimpleMCReconstruction.py and CostumHybridDetector.py scripts should work with the newest changes, they broke because the interpolater package was not installing properly.

However, some new tests with the new scripts you added would indeed be good. Though I would say this goes in a new PR.

…dDetector script

This reverts commit 1626f9e because it is in the unit tests.

…ad of interpolation. Remove import statements

lpyras · 2025-01-30T02:52:04Z

I tried to clean up. I am happy to approve this PR now.

MijnheerD

Why is it necessary to check for the number of observers and block interpolation with less than 50 antennas? If someone wants to try interpolation with few observers (eg for testing), they should be able to. I don't think it is up to this module to perform these checks, I think we can leave this responsibility to the user.

MijnheerD · 2025-01-30T09:48:19Z

NuRadioReco/modules/io/coreas/coreasInterpolator.py

@@ -303,6 +303,13 @@ def initialize_efield_interpolator(self, interp_lowfreq, interp_highfreq, **kwar

            return None

+        if self.obs_positions_showerplane[:, 0].shape[0] < 50:


Where does the number 50 come from?

I agree, that does not make much sense to me.

It’s not a great solution. I did that because I tried the module with a non-star shape pattern coreas file. So ideally one would check if the observer positions are a star shape pattern.

MijnheerD · 2025-01-30T09:50:02Z

NuRadioReco/examples/cr_analysis/sim_cr_single_station.py

 from NuRadioReco.utilities import units
 from NuRadioReco.detector.detector import Detector
 import NuRadioReco.modules.io.coreas.readCoREASStation
 import NuRadioReco.modules.efieldToVoltageConverter
 import NuRadioReco.modules.RNO_G.hardwareResponseIncorporator
 import NuRadioReco.modules.channelGenericNoiseAdder
-import NuRadioReco.modules.channelGalacticNoiseAdder
+#import NuRadioReco.modules.channelGalacticNoiseAdder


If the import statement is not necessary, please remove it entirely. The module exists though, so why remove this? Was it causing problems for you?

The galactic noise adder has a lot of dependencies. I thought for an example it is better if you don’t have to deal with them.

MijnheerD · 2025-01-30T09:52:09Z

NuRadioReco/modules/io/coreas/coreas.py

@@ -540,7 +540,7 @@ def add_electric_field_to_sim_station(

 def calculate_simulation_weights(positions, zenith, azimuth, site='summit', debug=False):
    """
-    Calculate weights according to the area that one simulated position in readCoreasStation represents.
+    Calculate weights according to the area that one observer position in a start shape patter represents.


MijnheerD · 2025-01-30T09:53:08Z

NuRadioReco/modules/io/coreas/coreas.py

+    if n_positions < 50:
+        logger.warning(
+            f"Only {n_positions} observer positions given to calculate weights. This is likely not a star shape pattern. Weights are set to 1.")
+        return np.ones(n_positions)


Similar to a previous comment, why is this should a deal breaker? And where does the number 50 come from?

The weight function is made for star shape pattern coreas files. It doesn’t really work for other observer positions. It’s not great solution I agree.

MijnheerD · 2025-01-30T09:57:01Z

NuRadioReco/modules/io/coreas/readCoREASStation.py

@@ -53,13 +54,14 @@ def run(self, detector):
            corsika_evt = coreas.read_CORSIKA7(input_file)
            coreas_sim_station = corsika_evt.get_station(0).get_sim_station()


Do we still need the sim station if the zenith/azimuth come from the shower now?

fschlueter · 2025-01-30T11:46:01Z

One comment: For the future, it might be useful to add a module which can export the interpolated simulation to a CoREAS style hdf5 file?

lpyras · 2025-01-30T15:44:06Z

Why is it necessary to check for the number of observers and block interpolation with less than 50 antennas? If someone wants to try interpolation with few observers (eg for testing), they should be able to. I don't think it is up to this module to perform these checks, I think we can leave this responsibility to the user.

Okay, I get your point. I stumbled over this, because I tried the examples with the example data, that lead to a lot of not understandable errors. Should we make a check for star shape patterns? Otherwise we could use a different example shower?

lpyras and others added 9 commits May 31, 2024 20:18

add docstrings

bc6eec4

refactor coreas interface

f7c9afb

add module to reconstruct efield from othorgonal antenna pairs

d1f6567

finish refactor, bugfix (group_id was used as channel_id to get anten…

cf1e900

…na position) code seems to give expected results now.

set efield position, fix efield reconstructor noise window definition

9a40073

Merge branch 'develop' into refactor_coreas_v3

c389b92

change corsika to corsika_evt to avoid confusion with open hdf5 file,…

88fc53f

… improve docstring

fix coreas.py

4931502

cg-laser and others added 7 commits June 6, 2024 12:44

adapt logging, fix fluence calculation

1af1f0d

allow to set declination

eefac87

remove tests and examples which rely on readCoreas. New examples and …

627ac86

…test will be part of PR #689

fix docstrings

1b041a8

Merge branch 'develop' into refactor_coreas_v3

33617b0

add automatic data download

58f7655

update readme

86f60de

lpyras added 2 commits June 18, 2024 17:04

example scripts for cr analysis with a detector array and a single st…

2768e16

…ation.

set correct start time for trace

b1578b1

lpyras added 4 commits August 5, 2024 17:57

write new functions which rely on readCorsika7()

524554c

refactor readCoREASShower to use readCOREAS7() function

42c3f59

rewrite readCoREASStation such that it uses readCorsika7()

acc8c43

delete old examples which is no longer applicable. Add new example code

63a8ddd

MijnheerD and others added 4 commits January 29, 2025 12:16

Change function name for clarity

5e8bf81

Use helper function

9ca6cf1

Update

0c27e10

Add deprecation warning

2126ffe

fschlueter previously approved these changes Jan 29, 2025

View reviewed changes

Readded the cr-pulse-interpolator dependency that went missing

d21e244

MijnheerD dismissed fschlueter’s stale review via d21e244 January 29, 2025 12:21

MijnheerD and others added 11 commits January 29, 2025 17:01

fixup! Use helper function

6ea6e95

use readCoREASStation instead of old module readCoREAS in CostumHybri…

f9cf4da

…dDetector script

add docstrings to cr examples

894a7df

minor docstring typos

8d5ea3f

rm outdated README in examples

a47d017

rm old air_shower examples since it doubles with cr_sim scripts

1626f9e

fix typo

1f20050

Revert "rm old air_shower examples since it doubles with cr_sim scripts"

9e85ea0

This reverts commit 1626f9e because it is in the unit tests.

minor changes to make readCoREASStation and the example work

12eddb6

improve det_array example

7ce6e87

catch hdf5 files with few observers and return closest observer inste…

7859d84

…ad of interpolation. Remove import statements

lpyras previously approved these changes Jan 30, 2025

View reviewed changes

MijnheerD reviewed Jan 30, 2025

View reviewed changes

Correct time units

1b4b7fc

MijnheerD dismissed lpyras’s stale review via 1b4b7fc January 30, 2025 10:10

squash! Correct time units

c8721c8

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Refactoring coreas to include interpolation #688

Refactoring coreas to include interpolation #688

cg-laser commented Jun 5, 2024

cg-laser commented Jun 5, 2024

cg-laser commented Jun 14, 2024

cg-laser commented Jun 14, 2024

lpyras commented Aug 5, 2024

lpyras commented Aug 5, 2024

lpyras commented Aug 6, 2024 •

edited

Loading

MijnheerD commented Aug 7, 2024

anelles commented Aug 7, 2024

fschlueter commented Jan 29, 2025

MijnheerD commented Jan 29, 2025

MijnheerD commented Jan 29, 2025 •

edited

Loading

lpyras commented Jan 30, 2025

MijnheerD left a comment •

edited

Loading

MijnheerD Jan 30, 2025

fschlueter Jan 30, 2025

lpyras Jan 30, 2025

MijnheerD Jan 30, 2025

lpyras Jan 30, 2025

MijnheerD Jan 30, 2025

MijnheerD Jan 30, 2025

lpyras Jan 30, 2025

MijnheerD Jan 30, 2025

fschlueter commented Jan 30, 2025

lpyras commented Jan 30, 2025

		@@ -303,6 +303,13 @@ def initialize_efield_interpolator(self, interp_lowfreq, interp_highfreq, **kwar

		return None

		if self.obs_positions_showerplane[:, 0].shape[0] < 50:

		@@ -53,13 +54,14 @@ def run(self, detector):
		corsika_evt = coreas.read_CORSIKA7(input_file)
		coreas_sim_station = corsika_evt.get_station(0).get_sim_station()

Refactoring coreas to include interpolation #688

Are you sure you want to change the base?

Refactoring coreas to include interpolation #688

Conversation

cg-laser commented Jun 5, 2024

cg-laser commented Jun 5, 2024

cg-laser commented Jun 14, 2024

cg-laser commented Jun 14, 2024

lpyras commented Aug 5, 2024

lpyras commented Aug 5, 2024

lpyras commented Aug 6, 2024 • edited Loading

MijnheerD commented Aug 7, 2024

anelles commented Aug 7, 2024

fschlueter commented Jan 29, 2025

MijnheerD commented Jan 29, 2025

MijnheerD commented Jan 29, 2025 • edited Loading

lpyras commented Jan 30, 2025

MijnheerD left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fschlueter commented Jan 30, 2025

lpyras commented Jan 30, 2025

lpyras commented Aug 6, 2024 •

edited

Loading

MijnheerD commented Jan 29, 2025 •

edited

Loading

MijnheerD left a comment •

edited

Loading