Starshapes clarity kit #65
Conversation
…rsika input system. Default stays the same as before
xrange to range
…radiotools into starshapes_clarity_KIT
… way of writing into list file
…radiotools into starshapes_clarity_KIT
There was a problem hiding this comment.
Thanks a lot for these improvements. I think it makes these scripts more usable, especially for newcomers.
I added a few inline comments that should be addressed before merging. I also saw that there is a conflict with the master branch. Can you merge the current master into this branch git merge master and resolve the conflicts?
A general comments about units. We don't enforce it but we are using the default units from NuRadio* https://github.com/nu-radio/NuRadioMC/blob/develop/NuRadioReco/utilities/units.py (which should be the same as the OFFLINE and GEANT default units).
Whenever no unit is specified, it should be in these default units. And ideally we always use these default units.
The only exception is "Grammage" which is defined in g/cm^2 due to historical reasons and because it gets unreasonable values when defined in the NuRadio default units.
| return corsika_azimuth | ||
|
|
||
|
|
||
| def get_rmax(X): |
There was a problem hiding this comment.
This will be very observation level and to a lesser extent also bandwidth dependent. Can you specify for which settings you estimated the relationship? To improve it for future use, the function could take the observation level as additional input which would allow to add additional parameterizations for other sites later.
There was a problem hiding this comment.
The input for the function is already observation level dependent since it uses the grammage at observation level when it's used in the next function . The grammage is calculated with "get_atmospherere" from the models.py file.
In any case, it needs more documentation.
I agree it's not great solution for automatically generating antenna rings. When using this branch to generate simulations recently, I also noticed that the new function "get_starshaped_pattern_radii" that I added has a dependency to a function from @fschlueter's RadioAnalysis framework to estimate the Cherenkov radius of a shower.
I used these two functions to, independly from geometry, generate antenna rings that go to big axis distances, but have a higher density from the shower core to a little further out than the Cherenkov radius.
I think it's an important functionality to have, but I might have to implement it differently. If I kept it mostly the same, then I need to implement the Cherenkov radius functions into radiotools as well.
…radiotools into starshapes_clarity_KIT
…radiotools into starshapes_clarity_KIT
…radiotools into starshapes_clarity_KIT
| from past.builtins import xrange | ||
| from radiotools.atmosphere import models | ||
| import sys | ||
| from miniradiotools.utils.cherenkov_radius import get_cherenkov_radius_model_from_depth |
There was a problem hiding this comment.
Is that a new external dependency
There was a problem hiding this comment.
yeah, I noticed it after the original changes and imported it from another another repository in order to use this branch short-term, but I need to fix it so it isn't using anything external. See my last message about the new functions
|
Everything should be addressed now. |
There was a problem hiding this comment.
I didn't go through all changes again but I trust that you have checked things properly, so let me approve it. Maybe @fschlueter want's to have another look?
| # old: cherenkov_angle_from_density_refractivity(rho, dxmax, n_asl, rho_0, ...) | ||
| # param of the cherenkov angle from star-shape simulations | ||
| # here cherenkov radius refers to radius of strongest emission | ||
| # is used to determine rmax in the RdIdealGrid simulations! | ||
| def cherenkov_angle_param( | ||
| height, dist, n0, model, | ||
| a=9.48990456e-01, b=4.48698860e+03, | ||
| c=1.43097665e+00, d=2.46630811e+06): | ||
|
|
||
| n = atm.get_n(height, n0=n0, model=model) | ||
| A = a - (b / dist) ** (c) - dist / d | ||
|
|
||
| return A * cherenkov_angle_model(n) | ||
| # # old param | ||
| # def cherenkov_angle_from_density(x, A=0.24905864, B=0.92165234): | ||
| # return np.deg2rad(A * np.log(x) + B) |
There was a problem hiding this comment.
Is that relevant to keep? I am not sure if that is worth keeping in here. You could also remove the "model" from all the function names?
There was a problem hiding this comment.
In my code I added model to the function names because I implemented the revent functions with the same name, I believe.
Anyway, I think it is a good idea to add them here.
There was a problem hiding this comment.
the "old" section is irrelevant and shouldn't be included. I missed it when copying the other functions.
I'll remove the model parts of the function names too, they seem unneccessary here.
| # calculate maximum distance of antenna from shower core (in shower plane) | ||
| # uses rough, hardcoded parametrisation | ||
| maxX = at.get_atmosphere(zenith, obs_level) | ||
| rmax = get_rmax(maxX) |
There was a problem hiding this comment.
Where is this coming form? Is that a (very) old parameterisation from Christian based on vertical showers? I would tend to define a maximal radius by a multiplicity of cherenkov radii.
There was a problem hiding this comment.
I think so. It's effectively doing that by generating antenna rings that stop about 5 Cherenkov radii out. I just did that with a little trial and error. I agree though that it's more elegant to use the Cherenkov radius directly, especially with the new functions I added.
I'll make the changes soon.
| obs_level=1400.0, | ||
| obs_level_corsika=None, | ||
| ground_plane=True, | ||
| Auger_CS=True, |
There was a problem hiding this comment.
That was one of the main reason I started working on this. The old code implicitly assumed that you are using the Auger coordinate system for giving the input parameters and just rotated them to the system CORSIKA uses somewhere in the function.
The main change I made is make this explicit and documented, and add the option to also use CORSIKA's coordinate system for the input. The Auger_CS option determines the input coordinate system.
There was a problem hiding this comment.
Sure, may comment was regarding all the capital letters in the variable name.
There was a problem hiding this comment.
good point, I'll make them lower case
Changes to models.py:
added 2 atmosphere models for experimental sites in China for GRAND
removed a deprecated module that didn't run for me (xrange -> range)
slightly widened a range in an atmospheric height calculation because it crashed for the 50-200 MHz frequency band
Changes to coordinatesystems.py:
just a typo in a print command
Changes to generate_coreas_sim.py:
extensively expanded the "write_list_star_pattern" function with detailed documentation. Everything new is at least briefly described. The default output of the function stays the same.
The main change is to replace the undocumented rotation by -90 degrees to transform to the Auger coordinate system with an explicit rotation matrix that is well documented. It also enables generating starshapes in Corsika's native coordinate system with a new input option: Auger_CS. By default it is set to True and performs the same -90 degree rotation as the current version of the code.
Error messages are included to catch inputs in the wrong units for angle and distances.
The function also now returns the azimuth angle that has to be written in the Corsika input to get the right shower direction for the corresponding starshape.
Added another optional input option to the function that writes a second .list file containing antenna positions in the shower plane coordinate system to enable visual checks of the starshape.
Added new functions to generate pre-made antenna rings with variable density that can conform to the shower geometry.