Refactoring coreas to include interpolation #688
Conversation
…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
|
@MijnheerD @lpyras I'm getting these warning from the interpolator module all the time: are we not using the interpolator correctly? |
…test will be part of PR #689
|
@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? |
|
The output of the reconstruction looks like this.
|
@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? |
I tested it before the coreas refactoring for effective area calculation. That yielded reasonable results. |
|
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. :) |
|
Sorry for coming in late on this, I thought I already answered some the questions but apparently I didn't.
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
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. |
|
I just tested the branch and it actually fails with a pip installed cr_interpolator on: |
| return obs_positions_vxB_vxvxB | ||
|
|
||
|
|
||
| def read_CORSIKA7(input_file, declination=None): |
There was a problem hiding this comment.
Can someone help me understand why the declination parameter is necessary here? Does it not come from the input_file itself?
There was a problem hiding this comment.
Without having checked the code and what the function does with the declination: No the declination does not come from the input_file. The normal corsika does not care about the declination at all as the coordinate system is defined by the magnetic north. There is the possibility to store the declination as meta data in an *reas file but it wont effect the simulation at all (its only intended to be stored along a simulation to use it once the simulation is converted to a different coordinate system in an reconstruction framework). Did that help?
There was a problem hiding this comment.
I see. But one could get it from the radiotools helper function, which has the 3D magnetic field, right? Do I understand that correctly?
I think Sjoerd added it at some point. But directly the link to the GitHub rather than the pip installation. |
… for important attributes
There was a problem hiding this comment.
I redid major parts of the code and fixed up (almost) all of the docstings, so type hints should now work better as well. The major things that remain are the following:
- The fluence interpolation is difficult to get working, because one has to set a fluence parameter manually. I think we should have an option for the fluences to be calculated from the traces (maybe one provides a function that does the operation, so that you can choose how to do this).
- I don't like the fact that parameters such as
zenith,azimuth, etc. are stored in multiple places (SimShower,SimStationand theElectricFields). I think we should make a choice for each of them, and always retrieve it from that place. - I have not checked the
calculate_simulation_weights()function, because I have no idea what it is supposed to do. - Should the plotting functions really be part of this module?
I also left some smaller comments in the code itself.
| sim_station.set_parameter(stnp.azimuth, azimuth) | ||
| sim_station.set_parameter(stnp.zenith, zenith) | ||
| sim_station.set_parameter(stnp.cr_energy, energy) | ||
| sim_station.set_magnetic_field_vector(magnetic_field_vector) |
There was a problem hiding this comment.
Personally I don't like the fact that we are storing this information twice (or maybe even three times), namely on the SimShower and SimStation level. Would it not be better if we just stored everything in the SimShower and use the SimStation only to store the traces?
There was a problem hiding this comment.
I have no strong opinion. I though with read_CORSIKA7() all comes from one place. I guess @cg-laser has to say something about this.
| "No high-level quantities in HDF5 file, not setting EM energy, this warning will be only printed once") | ||
| warning_printed_coreas_py = True | ||
| if coreas_sim_station.is_cosmic_ray(): | ||
| sim_station.set_is_cosmic_ray() |
There was a problem hiding this comment.
Should we keep this? This will always be True if the Event is created with read_CORSIKA7
There was a problem hiding this comment.
It's also legacy. @cg-laser ?
| efield_pos = [] | ||
| for efield in efields: | ||
| efield_pos.append(efield.get_position()) | ||
| efield_pos = np.array(efield_pos) |
There was a problem hiding this comment.
This is an array shaped like (nr_of_efields, 3), so it makes no sense to use this for the core. But I don't understand what the core is supposed to be here? Just the absolute position of the station? Would it not be better to remove the detector and station_id parameters and just allow the user to specify the core they want?
There was a problem hiding this comment.
This is from before we had the interpolation. The idea was, that one uses a shower n_obs times to calculate the effective area. So to avoid interpolation we place the core in a way, that a certain observer is on top of a single detector station. I agree, that this is confusing but it serves its purpose. For this, also the calculate_simulation_weights() function was used. I will follow up in another comment.
|
|
||
| return self.fluence_interpolator | ||
|
|
||
| def plot_footprint_fluence(self, radius=300): |
There was a problem hiding this comment.
Should this be part of the interpolation module?
There was a problem hiding this comment.
I found it useful, but I don't insist on it :)
| self.electric_field_on_sky, | ||
| signals_start_times=self.efield_times[:, 0], | ||
| lowfreq=(interp_lowfreq - 0.01) / units.MHz, | ||
| highfreq=(interp_highfreq + 0.01) / units.MHz, |
There was a problem hiding this comment.
I added a margin to avoid edge effects when applying filters. But I am not sure if it is really necessary. Feel free to delete it.
| """ | ||
| fluence_per_position = [ | ||
| np.sum(efield[quantity]) for efield in self.sim_station.get_electric_fields() | ||
| ] # the fluence is calculated per polarization, so we need to sum them up |
There was a problem hiding this comment.
There should be an option that this is calculated from the traces?
There was a problem hiding this comment.
@cg-laser ? I am happy with everything.
| else: | ||
| return self.geo_star_radius | ||
|
|
||
| def initialize_efield_interpolator(self, interp_lowfreq, interp_highfreq): |
There was a problem hiding this comment.
With the LOFAR version we had an optional parameter for the interpolation parameters. I think that was a clean solution to have some flexibility for the user. What do you think?
There was a problem hiding this comment.
Can you post a link to the LOFAR interpolation?
| return fig, ax | ||
|
|
||
| def get_interp_efield_value(self, position_on_ground, core): | ||
| """ |
There was a problem hiding this comment.
I am removing the core parameter, because I think it is very confusing. This module is purely about the interpolation,so everything should be from the reference frame of the original starshape IMHO. If one wants to shift the antenna positions, to account for a core shift, I think the should do this higher up and subtract the core from the position_on_ground value. What do you think?
There was a problem hiding this comment.
sure, we can move the function to another place. I placed it here because the interpolator object is stored in that module. Ideally, it would live in readCoREASDetector.py I guess.
|
I just added one important TODO item to the |
|
Thank you Mitja for putting in all the work :)
We used this function for the effective area calculation. Because of limited resources we only used one detector station and used the simulated showers multiple times. To avoid interpolation, we placed the shower core in a way, that the detector station is at the same position as the observer and determined the trigger efficiency. To figure out which area of the shower we would detect with that station we calculated a "cell" around that station that is represented by that observer position. For the effective area we summed up the cells on which we triggered for each energy and zenith bin. The function 
For a whole detector array we used a different method, also to simulated coincidences. But as a first guess this works fine. I think we should keep it, since we are also keeping the readCoREASStation.py. |
| @@ -267,41 +267,7 @@ def read_CORSIKA7(input_file, declination=None): | |||
| zenith, azimuth, magnetic_field_vector = get_angles(corsika, declination) | |||
| energy = corsika['inputs'].attrs["ERANGE"][0] * units.GeV # Assume fixed energy | |||
|
|
|||
| # Create RadioShower to store simulation parameters in Event | |||
There was a problem hiding this comment.
I am not sure if @cg-laser goes along with that. The idea was, that read_CORSIKA7() is the "mother" function from which all other functions are derived and which can easily be replaced by read_CORSIKA8()
| @@ -426,63 +392,62 @@ def create_sim_station(station_id, evt, weight=None): | |||
| return sim_station | |||
|
|
|||
|
|
|||
| def make_sim_shower(corsika, observer=None, detector=None, station_id=None): | |||
| def make_sim_shower(corsika, declination=0): | |||
There was a problem hiding this comment.
I think this function should have another name to allow for a deprecation warning. As mentioned before, I think it should be part of read_Corsika.
| shower_axis_vector = hp.spherical_to_cartesian(zenith / units.rad, azimuth / units.rad) | ||
| geomagnetic_angle = hp.get_angle(magnetic_field_vector, shower_axis_vector) * units.rad |
There was a problem hiding this comment.
I understand you want to be "super safe" here. But I think the idea of the unit system of NuRadio is that once we read external data in we do not have to use units ever again until we start plotting or exporting data again. This also means we will never change the default units.
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.