Atmosphere models for h_max to X_max conversion #2000
Conversation
|
Check out this pull request on See visual diffs & provide feedback on Jupyter Notebooks. Powered by ReviewNB |
Codecov ReportBase: 92.50% // Head: 92.51% // Increases project coverage by
Additional details and impacted files@@ Coverage Diff @@
## master #2000 +/- ##
==========================================
+ Coverage 92.50% 92.51% +0.01%
==========================================
Files 197 199 +2
Lines 16311 16561 +250
==========================================
+ Hits 15089 15322 +233
- Misses 1222 1239 +17
Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here. ☔ View full report at Codecov. |
kosack
force-pushed
the
feature/atmosphere
branch
from
232e8be to
14ad684
Compare
|
Hi Max,
I am not aware of simtel files from a single run with more than
an I/O block of type 1216 - if you found such a case, let me
know how it was produced. If they are identical, using the
first one is OK but there are situations possible where a
first one might just contain part of the information and
might later get complemented by a second such block
with the other part. Typical CORSIKA simulation should
fill both parts right away though.
So, well, there are two parts of content in this data block.
One part, which you probably want to convert, is the copy
of the atmprof*.dat file read from CORSIKA. The other part
is the result of the 5-layer fit used by the CORSIKA EGS part.
The code running with CORSIKA should report both parts,
if available, but you might end up with just one part.
The table part would be missing if you run corsika without
an atmospheric profile table. The other ways to configure
the atmospheric profiles in CORSIKA only set the 5-layer
parameters. In that case there is no table to report.
The thickc() and heighc() functions in mc_atmprof.c are
C implementations of what CORSIKA is doing internally
with the 5-layer atmosphere models.
Too bad that EGS part is hard-coded to use the 5-layer
function but it cannot accurately describe the actual
atmosphere. If you want to lookup some height from an
atmospheric depth in the longitudinal profiles, you first
need to make up your mind which description is the more
appropriate one for the task. Since most Cherenkov light
gets emitted from e+-, the 5-layer description might actually
be closer to what CORSIKA originally generated. Xmax is a
mixed beast, at least for hadron showers, so there would be
no really correct way. Hadrons and muons follow directly
the table profile, if any got set up and you don't run a
CURVED option CORSIKA.
Note that interpolation of the tabulated profiles is quite non-trivial.
In particular, you want to have height(thick(h)) == h, at an
acceptable accuracy.
Linear interpolation in a mostly exponential density profile
will obviously give wrong answers, even though that looks
reversible. The atmo.c code goes through some effort to
cspline interpolate in log(thick) versus height - except at the
top where thick=0 is not suitable for the log while csplines
are not very forgiving for any sudden transitions (jumps in
second derivative). The reversion through the heighx_ and
thickx_ functions in atmo.c are good to about 0.2 mm throughout
the whole atmosphere while direct interpolation methods
might easily miss by a 100 meters.
Best regards,
Konrad
Am 11.07.22 um 18:08 schrieb Maximilian Nöthe:
… @maxnoe commented on this pull request.
> +
+
+def read_simtel_profile(simtelfile: str) -> Table:
+ """
+ Read an atmosphere profile from a SimTelArray file as an astropy Table
+
+ Returns
+ -------
+ Table:
+ table with columns `height`, `density`, and `column_density`
+ along with associated metadata
+ """
+ import eventio
+
+ tables = []
+ with eventio.SimTelFile(simtelfile) as simtel:
@kbernloehr Any reason why there could be more than one atmospheric profile in a simtel file?
|
|
Thanks. As it's not intended for simulation, the use case for this package probably doesn't require extreme precision in the atmosphere profile or it's interpolation as other systematics dominate, but in any case I'll switch to a log interpolation and use cubic splines. I'll also add a disclaimer in the documentation! |
|
Hi Karl,
I agree that for Xmax lookup you can be a bit more
relaxed. You definitely want to leave out the last point
(thick=0) but perhaps even the second to last point
to avoid wiggles in the csplines. The shower maximum
will never be at these altitudes.
And, even though all of our available productions used
tabulated atmospheres, you might want to be prepared
for those where the atmospheric profile is defined by other
means (switch to 5-layer description).
Best,
Konrad
Am 12.07.22 um 10:11 schrieb Karl Kosack:
… Thanks. As it's not intended for simulation, the use case for this package probably doesn't require extreme precision in the atmosphere profile or it's interpolation as other systematics dominate, but in any case I'll switch to a log interpolation and use splines. I'll also add a disclaimer in the documentation!
|
kosack
force-pushed
the
feature/atmosphere
branch
from
b066b92 to
3cb4d29
Compare
Do you want to add a class for that description as well, @kosack ? As a fallback if the tabulated version is not available? |
|
Interpolation in the wiggly areas looks like this now (with cubic splines used):
|
Is that documented somewhere and do we have a file that has that format atmosphere? If so, I can add it, or we could always add a subclass in a later version |
|
The 5-layer coefficients are, of course, documented in the CORSIKA
User's Guide.
See, in particular, Appendix F Atmospheres. The 'DATM' variable is not
documented
there and is just 1/CATM, unless CATM=0. Used in CORSIKA because
multiplication
is faster than division.
Am 13.07.22 um 15:43 schrieb Karl Kosack:
…>> means (switch to 5-layer description).
> Do you want to add a class for that description as well, @kosack ? As a fallback if the tabulated version is not available?
Is that documented somewhere and do we have a file that has that format atmosphere? If so, I can add it, or we could always add a subclass in a later version
|
All the files should have it, the data array is next to the tabulated values in |
|
I added a
|
|
The last item to implement here is the i/o. Since there could in principle be more than one version of the atmosphere profile in the file (one from the simulation, and one from real measurements for observation data), I think the ones loaded by the SimTelEventSource will go in /simulation/service/atmosphere_profile, which then allows one to add a /dl2/service/atmosphere_profile later for the frequently updated/measured one. |
I don't understand why this split is necessary. Simulated files only have the simulated atmosphere, observed files only have a measured one. Why can we not use the same location? And wouldn't the atmosphere profile be monitoring not service information? |
There seems to be a discontinuity in the first plot, reading the CORSIKA manual, this should not happen? |
|
Hi Max, A jump in density looks ugly but could even kind of happen in a The 5-layer thickness should never be far from the tabulated Since the jump looked too big for my gut feeling, I look into Atmospheric profile 99 to be read from file 'atmprof_ecmwf_south_winter_fixed.dat'. Atmospheric profile 99 with 55 levels read from file atmprof_ecmwf_south_winter_fixed.dat Initialize atmosphere ranging from 0.000 to 120.000 km a.s.l. Results of the atmosphere fit: Altitude [km] rho(table) rho(fit) thick(table) thick(fit) The border boundaries are identical to what your earlier post showed, |
I only had the bare array from the simtel file, so I just guessed at what it meant (following the Corsika manual), however I think it looks ok, but I could have a index-1 bug somewhere. The thickness profile looks fine (very little difference from a table profile), and the jumps in the density are just because I simply used the same parameters for thickness but took their derivative to get density, so those are just cases where the derivative wasn't continuous between bins I guess (or perhaps a bug - i'll check). I interpreted the data as a table as follows: |
|
The comparison looks like this: |
It's just a bookkeeping thing... one might want to eventually have both the simulated and measured profiles in the same file, for comparison reasons, so having them in different places of the data model could help that. But for now, I'll just concentrate on the simulated version, so I'll add it to |
|
It looks like my implementation is identical to @kbernloehr : fit_konrad = np.array(
[
[0.00*100000, -140.508, 1178.05, 994186,0],
[9.75*100000, -18.4377, 1265.08, 708915,0],
[19.00*100000, 0.217565, 1349.22, 636143,0],
[46.00*100000, -0.000201796, 703.745, 721128,0],
[106.00*100000, 0.000763128, 1, 1.57247e10,0],
]
)
k5 = FiveLayerAtmosphereDensityProfile.from_array(fit_konrad)
So I thikn the implementation is ok. [Edit: I had incorrect binning in h in my previous version of this plot] |
replaced by ABC
…nto feature/atmosphere
and change supported table versions to a set
| @@ -1,3 +1,5 @@ | |||
| """ tests of SimTelEventSource """ | |||
| # pylint: disable=import-outside-toplevel | |||
There was a problem hiding this comment.
For another PR: can we disable import-outside-toplevel for all test files?
There was a problem hiding this comment.
I haven't found a way to filter by filename, but it would be nice! We could just disable it globally since it's used sparingly and correctly for the most part.
Features
Adds an
AtmosphereDensityProfileclass hierarchy with three implementations:ExponentialAtmosphereDensityProfile(purely mathematical)TableAtmsophereDensityProfile(cubic interpolation of a tabulated model)FiveLayerAtmosphereDensityProfile(implements the Corsika 5-layer representation where each layer is either an exponential or linear model)This adds a new
atmosphere_density_profileattribute to EventSource, with an implementation inSimTelEventSource(reading them from simtel files) andHDF5EventSource, as well as inDataWriter(propagating them to output files).By default SimTelEventSource will return a TableAtmosphereDensityProfile if it is found in the data, otherwise it falls back to a FiveLayerAtmosphereDensityProfile, or None if no profile exists. This is user-selectable with the
SimTelEventSource.atmosphere_model_choiceoption (AUTO, TABLE, FIVELAYER)Example:
Loading one from a simtel or ctapipe file can be done with an EventSource:
All profiles can also be plotted: