search

cta-observatory / ctapipe

Low-level data processing pipeline software for CTAO or similar arrays of Imaging Atmospheric Cherenkov Telescopes
https://ctapipe.readthedocs.org
BSD 3-Clause "New" or "Revised" License
62 stars 266 forks source link

Is the line_of_sight_integral function really correct? #2515

Closed maxnoe closed 3 months ago

maxnoe commented 4 months ago 
          Another question about the functions currently there: Is the line_of_sight_integral function really correct? Because the distance you put into it is a distance from the observer, but the integral functions uses heights a.s.l.

Originally posted by @gschwefer in https://github.com/cta-observatory/ctapipe/issues/2422#issuecomment-1790265905

kosack commented 4 months ago

Sorry, didn't read the full question correctly. Yes, we should check that distance to observer is used correctly. I corrected it in the h_max code, I think, but likely there is something still missing (like the distance from the observer projected along the LOS to sea level)

kosack commented 4 months ago

It think the best is just to change the help string to be correct: it is not the distance from the observer, but from sea level. It's actually the height above sea level of the point you want to turn into a grammage. Then it's consistent witth how we now define h_max (which is from sea level). From that height, you integrate out to infinity along the zenith direction. Calling it "height_to_slant_depth" would probably be more clear as well. I'll try to do a full cleanup of this shortly, and also integrate #2422 . See also #2403 (merged), which fixed the definiion of h_max, which was originally relative to the obervatory, but is not correctly from sea level

gschwefer commented 4 months ago

I'm having a look now, also at #2422

gschwefer commented 4 months ago

Alright I looked at this now...The circle (vertical) overburden ->height a.s.l. ->(vertical) overburden is fully implemented via the functions "integral" and "height_from_overburden". These are what's actually used in all the internal calculations. The functions "height_from_slant_depth" and "line_of_sight_integral" are only convenience functions to account for the zenith angle automatically and not used in any internal calculations. I agree with you, Karl: I will rewrite the "line_of_sight_integral" so it's the proper inverse of "height_from_slant_depth" and integrate this in #2422

kosack commented 4 months ago

I think you should also refactor ImPACT to use the "convenience" functions, since there is an implicit assumption on the atmosphere shape that is used, and we may improve that in the future. So really, these are not convenience functions: they implement a specific atmosphere model that may change in a future improvement, and thus all Reconstructors should use those functions directly and not implement their own. For example, we know currently for zenith angles > 70°, this model is incorrect. CORSIKA employs a trick for these higher zenith angles, keeping a flat atmosphere, and we could implement that in the future, or even move to a fully spherical atmosphere, but doing so would break ImPACT if it has it's own internal implementation of how to compute slant depth.

gschwefer commented 4 months ago

Yes, I will of course address these things in ImPACT. So if you want to keep this issue open for that, that's alright, but it should not block the merging of #2422, right?

kosack commented 4 months ago

Yes, I will of course address these things in ImPACT. So if you want to keep this issue open for that, that's alright, but it should not block the merging of https://github.com/cta-observatory/ctapipe/pull/2422, right?

No, that can be a separate PR of course! It was just to point out that reconstructors should not implement their own conversion from slant to height and just use the ones from the atmosphere module.

By the way, reposting this image from the previous issue/discussion, since it illustrates this issue, where the distance to the maximum should be along the shower-axis from sea level. image

Just note that the "line of sight integral" (x'(max) here) is from the shower-max point to infinity, so he distance from the observer doesn't really matter (as long as you properly compute the shower-max point. In the current implementation (after #2403), it is correctly computed from sea level, so in principle everything should work ok. It is only the inverse-transform where you have to be careful.