search

nu-radio / radiotools

A tool package for cosmic-ray and neutrino radio detectors
GNU General Public License v3.0
8 stars 5 forks source link

Fix model calculation include obs lvl #45

Closed fschlueter closed 3 years ago

fschlueter commented 3 years ago

Hi, this branch is not jet ready to be merged but needs some discussion. This PR picks up #38 and issue #41.

Briefly: Functions like get_vertical_height, get_distance_xmax, etc which take a zenith angle as input calculate the height, grammage, etc for a point in space. This point in space in only properly defined with a observation level which defines the inclination, i.e., zenith angle.

I ran some cross-checks on this. Compare to the result provided bei CORSIKA/CoREAS. The improvement is not significant (it was already good before).

The plots below show shower with zenith angle between 65 and 85 binned in 2.5deg steps.

With the "bug" dmax_comparison_rt_corsika_xmax_bug

With "bug" fixed for numerical calculation (applied for zenith angle > 80 deg, default) dmax_comparison_rt_corsika_xmax_bug_fix

With "bug" fixed for numerical calculation (applied for all zenith angle) dmax_comparison_rt_corsika_xmax_bug_fix_num

fschlueter commented 3 years ago

The "taylor" calculation can not be "fixed" by just adding the obs_lvl to the function call. The problem here is that the expansion of the "a" parameter depending of the zenith angle is performed with the obs_lvl = 0. I think there are two ways to fix it:

  1. Repeat the expansion for a defined obs_lvl. This has the caveat that the model would depend on the obs_lvl. For each different obs_lvl one would need a new instance of the model. However it would a bit "safer" as the obs_level has only to be set once by the user.
  2. To get "a" from the cached value convert the zenith angle (at obs_lvl) to the zenith angle at sea level. Something similar is done in refractivity.py. One would need to ensure that this conversion is robust.
fschlueter commented 3 years ago

With the last commit also the calculation via taylor looks fine:

Use taylor for zenith angle < 80deg dmax_comparison_rt_corsika_xmax_bug_fix_taylor_too

cg-laser commented 3 years ago

Hey Felix, thanks for doing this work! Funnily today I actually noticed that the observation level is not taken into account in the get_atomosphere function. I briefly scanned your code and it looks good, but I mainly looked at your plots. They are pretty convincing. I agree that the best way is to just correct for the zenith angle on the fly. Can you make another plot of how big this correction is? E.g. for the auger observation level? And then zenith correction vs. zenith angle?

fschlueter commented 3 years ago

Hi @cg-laser, not quite sure which correction you want to see (alongside the effect of the zenith angle correction). In this plot the difference of get_vertical_height (for an depth of 700 g/cm2) with and without obs_level passed is shown along with the effect on the geometrical distance to this point and the difference in the zenith angle at 0m and 1400m. show_obs_lvl_bug

fschlueter commented 3 years ago

Hi @cg-laser, not quite sure which correction you want to see (alongside the effect of the zenith angle correction). In this plot the difference of get_vertical_height (for an depth of 700 g/cm2) with and without obs_level passed is shown along with the effect on the geometrical distance to this point and the difference in the zenith angle at 0m and 1400m. show_obs_lvl_bug

cg-laser commented 3 years ago

@fschlueter thanks, this is what I was looking for. So the zenith correction is at most 0.4deg at close to 90deg. This gives me a feel for the effect. I think the PR is ready to be merged?