h_max is inconsistent between HillasIntersection and HillasReconstructor

[gschwefer]

Both HillasReconstructor and HillasIntersection calculate a quantity h_max and return it in the ReconstructedGeometryContainer, but they are inconistent in the physical meaning of h_max: For HillasIntersection, it is the vertical height of the shower maximum above the array, for HillasReconstructor it is the distance of the shower maximum location from the array center, i.e. kind of a "slant height"

Whatever h_max really should be, it should of course be the same everywhere.

[gschwefer]

See the related discussion in #2305 about what quantity returned from ImPACT

[gschwefer]

Also, should the Hillas reconstructors in addition to h_max return the atmospheric slant depth xmax?

[kosack]

Yes, I believe we should return both h_max (which should really be the height of shower max, the z coordinate), and the distance to shower max (along the line of sight). Slant depth (what we call in the atmosphere module "column density") is then the distance transformed to a grammage using an atmosphere model, whereas x_max is the h_max transformed to grammage units.

And to be more explicit, we might want to rename these to e.g.:

• shower_max_distance (m)
• h_max → shower_max_height (m)

or even just shower_distance, shower_height, since we usually define the shower's centroid point as the "Shower maximium", i.e. where most of the cherenkov light is emitted

[gschwefer]

Ok, I have one question: Do I understand you correctly that you want to give the name x_max to the vertical grammage equivalent to shower_max_height ?

[maxnoe]

I think that would be bad, since e.g. Auger uses x max for the depth of the shower maximum and I'm pretty sure that's generally the more common usage.

[kosack]

We might want to avoid x_max altogther if that is the case. Corsika/SimTelArray define it as a height (mass overburden), not corrected for the zenith angle. If Auger uses a different definition, it's clear that it is not well defined and we should pick names that are obvious.

Also "depth" as a word usually does not contain the cos factor (think water depth in the ocean, which is not a line-of-site measurement), so it's not even clear to me that Auger's definition is different. Normally, "depth", "height" and "altitude" (when not an angle) all mean really just the z-coordinate of a point in a fluid/gas defined from different 0 points.

The path-length from the observer to a point in a medium is different: it is the slant-depth, column-density, etc. It implies an integration over a path, which in our case is a line.

[kosack]

What we are describing is a point - the shower maximum. Here is a summary:

• solid red line: integrating along a path that starts at infinity and ends at the shower maximum point is what we currently call the "shower depth" , sometimes called "mass overburden"
• solid black line: in distance units from the ground, this is the altitude or height of shower maximum
• dotted black line: the distance from the observatory to the shower max forms a line, the "line-of-sight", the distance from the ground to the shower max along this line is the "shower distance" (unless you have a better name?)
• dotted red line: integral along the line-of-sight from infinity to the shower maximum is the slant depth or column density

[maxnoe]

Corsika/SimTelArray define it as a height (mass overburden), not corrected for the zenith angle.

Actually it depends... CORSIKA has a compile time option for it (SLANT), if set, longittudinal development is sampled in slant depth, if not in vertical depth, which what determines the value of the true xmax, since that comes from a fit to that distribution.

I guess AUGER is of course always using SLANT. We are probably not, but maybe we should? At least for high zenith angles?

[kosack]

Yes, I mentioned that in the other issue: in CTA, we do not use the SLANT option. Probably it's too late to enable it.

In any case, we should just come up with clearer names than x_max... or at least be consistent. Is the SLANT option written to the simtel output in a way we can read it?

[kosack]

In the atmosphere module, the density profiles are all defined in terms of height from ground, and each defines an integral along the height direction (so everything is 1D), returning depth. The profile.line_of_sight_integral() does the 2D integration from infinity to a point, along the line-of-sight. It currently assumes a cartesian atmosphere for simplicity, but a proper numeric integral could be done for a spherical atmosphere or even a more realistic oblate spheroid if we really wanted, though the difference is small below 60° zenith angle. That's just an implementation detail, and probably not worth worrying about yet.

[kosack]

Note that neither of these are really what you want for doing particle physics: you want the integral along the line of the shower direction from infinity to the max point, not the line-of-sight from the observer. So even if we switch the definition, we still have the wrong answer. I mean the atmosphere module is still correct, you just have to give it the zenith angle of the shower not the observation

Though it's not really a problem: HillasReconstructor reconstructs the distance from the array to the shower max. But it also computes the full shower axis since it gets the impact point and the point-of-origin, so you can compute the proper slant depth from that. It's why we don't write out x_max currently... But h_max is probably the wrong name for that attribute, and we could simply just properly compute it as a height and be consistent. Then all the complexity of computing the correct x_max still has to happen externally

[maxnoe]

Also related: https://github.com/cta-observatory/ctapipe/issues/1030

[kosack]

Yes, I thought of tat as well. Here's a better diagram of what we have now...:

But really we should add the observatory height, since the atmosphere model is from sea level.

So we want this:

[kosack]

Anyhow, I suggest the following:

For the definition of height, there are two options:

1. define heights as relative to the observatory's center point (0,0,0), e.g. consistent with the current definition of the telescope heights, which are also relative to the array center point (in x,y,z), where a telescope at height 10m is 10m above the array plane. Then in the conversion to x_max, the suer has to be careful to add the observatory height.
2. add source.subarray.reference_location.geodetic.height to h_max, so heights are real heights above sea level. however, we could leave telescope z coordinates relative, which still makes sense and won't break the reconstruction.

Probably option 2 is safest.

Then fix the following bugs / add features:

• HillasReconstructor currently returns h_max as the distance to the shower_max from the array center (0,0,0). This is just wrong. we should just take the z-coordinate (and add the reference height) and then it is a true "height". If we want the distance in meters to be stored additionally, we could compute the distance from the impact point on the ground projected to sea level, but I'm not sure that is very useful for anything.
• define x_max as just the vertical depth of the shower max (could rename it to shower_depth as mentioned above to avoid confusion, at the expense of changing the data model)
• define the slant depth as the line integral from infinity to the shower max point along the line defined by the shower zenith angle (not the pointing one). Name that shower_slant_depth everywhere.
• make a helper function that uses atmosphere module and the shower reconstruction to compute the real shower depth and slant depth (and vice versa) so we have one place where it is done correctly
• Fix HillasIntersection which does a few bad things: it does add the observatory depth (but as a hard-coded number which is wrong!), and it also calls the depth in m "x_max" internally, which is even more confusing. It never actually does the conversion

If that sounds ok, I'll make a preliminary PR to fix at least this in the HillasReconstructor and HillasIntersection and update naming in the atmosphere module. So far nothing other than ImPACT actually uses x_max or h_max, so that can be fixed separately.

[maxnoe]

Sounds good to me

[kosack]

i did run into the issue that not all of our test files provide a subarray reference_location... So adding the height fails for those. Either we should update the example files, or have a fallback.

In particular:

lst_prod3_calibration_and_mcphotons.simtel.zst

[maxnoe]

i did run into the issue that not all of our test files provide a subarray reference_location... So adding the height fails for those. Either we should update the example files, or have a fallback.

The obs level height should be available for all and could be used as fallback. Specifying the reference location via simtel metadata is fairly new (Prod6 prototype files).

We could maybe fill the reference location with dummy lat/lon of 0 but the correct altitude from obslevel in that case: https://en.wikipedia.org/wiki/Null_Island

[kosack]

We could maybe fill the reference location with dummy lat/lon of 0 but the correct altitude from obslevel in that case: https://en.wikipedia.org/wiki/Null_Island

That sounds reasonable. I'll open a separate issue or PR for that. For now I fixed the problem by just switching what file is used for generating the dummy subarray in the tests.

[kosack]

@maxnoe Where can I access obslev? Is it only in corsika_inputcards, or is it parsed out somewhere?