Telescope z positions are raised 10.5 m higher

YoshikiOhtani commented 2 years ago

Just in case let me report one thing that I have found (sorry if it is already reported somewhere...).

In the CORSIKA input cards, the telescope positions are defined as follows:

https://github.com/cta-observatory/lst-sim-config/blob/793c54962b21841a19e0a84794f1b79f12d2e795/config_LST1_2MAGIC_diffuse_gammas.cfg#L83-L85

However, when I opened one of LST-PROD2 simtel files I get the positions whose z coordinates are a bit raised:

In [10]: source.subarray.tel_coords
Out[10]: 
<SkyCoord (GroundFrame): (x, y, z) in m
    [( -6.33599997,  60.40499878, 12.5       ),
     ( 41.05400085, -79.27500153, 10.75      ),
     (-29.45600128, -31.29500008, 11.92000008)]>

And then I looked into the log file and see the following sentences:

Number of simulated telescopes: 3 All telescopes are raised by 10.50 m Telescope 1 at x= -6.34 m, y= 60.41 m, z= 12.50 m with r=12.50 m Telescope 2 at x= 41.05 m, y=-79.28 m, z= 10.75 m with r=10.00 m Telescope 3 at x=-29.46 m, y=-31.30 m, z= 11.92 m with r=10.00 m

It seems that the positions are raised so that the entire spheres defined at the telescope positions are above z = 0. I think it doesn't harm anything, but it would be good to modify the positions in the input card to make them consistent between input and output.

moralejo commented 2 years ago

I don't understand, the coordinates here:

https://github.com/cta-observatory/lst-sim-config/blob/793c54962b21841a19e0a84794f1b79f12d2e795/config_LST1_2MAGIC_diffuse_gammas.cfg#L83-L85

are intended to be the centers of rotation of the mirror dishes (all within 1 m of the corresponding dish center). The fourth column is the sphere radius, which should completely contain the telescope. See https://forge.in2p3.fr/issues/35807#note-98

I don't think there is any need to shift the sphere to be above ground. Anyone knows why is that done?

maxnoe commented 2 years ago

I don't think there is any need to shift the sphere to be above ground. Anyone knows why is that done?

I think this is needed for the ray-tracing code going back from position on the ground to the position where the photon might have it the telescope, but not sure. I remember reading about this requirement before but couldn't find it when I now quickly searched in the corsika manual and the iact/simtel docs

jsitarek commented 2 years ago

Hi @moralejo. The code in corsika/bernlohr/iact.c checks if ztel - rtel < 0 and if so raises all the telescopes by the difference. The effect is as @YoshikiOhtani said that the entire sphere around the telescope is above the ground level. I guess this is to avoid a situation when e.g. a photon crosses the horizontal plane containing the telescope center at a distance of a bit more than the radius of the sphere (so the photon is lost and not tracked more), but the continuation of the path "underground" would cross the sphere so the photon should have been kept.

moralejo commented 2 years ago

Hi @moralejo. The code in corsika/bernlohr/iact.c checks if ztel - rtel < 0 and if so raises all the telescopes by the difference. The effect is as @YoshikiOhtani said that the entire sphere around the telescope is above the ground level. I guess this is to avoid a situation when e.g. a photon crosses the horizontal plane containing the telescope center at a distance of a bit more than the radius of the sphere (so the photon is lost and not tracked more), but the continuation of the path "underground" would cross the sphere so the photon should have been kept.

So you mean the photons are tracked only to ground level, and if the track has not intersected the sphere by then they are ignored. I don't understand why that requires to modify explicitly the input telescope coordinates, but ok... In any case this will make no difference whatsoever. I would be fine with eventually doing the change proposed by @YoshikiOhtani, but not now that we are in the middle of the all-sky production.

maxnoe commented 1 year ago

Investigating the transformation between lon / lat height values (i.e. an astropy.coordinates.EarthLocation) and the simulation coordinates (i.e. ctapipe.coordinates.GroundFrame) and discussing with @moralejo we concluded:

the issue here is a real issue that should be fixed
The telescopes are "actually moved" so the positions are now incorrect in the simulations
However, the effect of this on the actual simulated events will be negligible
We'd still like to fix this for the next round of productions, so that we can have the same and correct locations of the telescopes in both simulations and observed data

Due to how sim_telarray and CORSIKA work and the locations / properties of the telescopes I would thus propose the following:

Take the measured locations of the central pins of LST-1, MAGIC-1, MAGIC-2 from @moralejo's presentation here: https://redmine.cta-observatory.org/attachments/65827
Add the distance between the central pin and the axis of rotation (from blueprints, also in the same presentation):
- LST: 15.883 m
- MAGIC: 10.6 m

This gives us the locations in lat/lon/height of the telescope center of rotation:

name,lon,lat,height
LST-1,-17.89149701,28.7615261,2199.883
MAGIC-1,-17.8900665372,28.7619439944,2195.363
MAGIC-2,-17.8905595556,28.7613094808,2197.364

The observation level has to be lower than these points minus the radius of the light sphere of each telescope. So I would propose to put the observation level at:

OBSLEVEL / ALTITUDE 2180 m

To transform from the simulation frame to/from earth locations, I opened this PR in ctapipe: https://github.com/cta-observatory/ctapipe/pull/2167 and added the same code temporarily to ctapipe_io_lst here: https://github.com/cta-observatory/ctapipe_io_lst/pull/161

Here is the full script performing the calculations, both with astropy and with pyproj, which should be the same as the transformations done before using the cs2cs command line tool (which comes from proj and pyproj offers the python bindings for that):


import numpy as np

from astropy.coordinates import EarthLocation, AltAz, ITRS
import astropy.units as u
from pyproj import CRS, Transformer

from ctapipe_io_lst.ground_frame import ground_frame_from_earth_location

tels = [
    dict(
        name="LST-1",
        loc=EarthLocation(lon=-17.89149701 * u.deg, lat=28.76152611 * u.deg, height=2199.883 * u.m),
        area=23**2 * u.m**2,
    ),
    dict(
        name="MAGIC-1",
        loc=EarthLocation(lon=-17.8900665372 * u.deg, lat=28.7619439944 * u.deg, height=2195.363 * u.m),
        area=17**2 * u.m**2,
    ),
    dict(
        name="MAGIC-2",
        loc=EarthLocation(lon=-17.8905595556 * u.deg, lat=28.7613094808 * u.deg, height=2197.364 * u.m),
        area=17**2 * u.m**2,
    ),
]

areas = [tel["area"].to_value(u.m**2) for tel in tels]
ref_lon = np.average([tel['loc'].lon.deg for tel in tels], weights=areas)
ref_lat = np.average([tel['loc'].lat.deg for tel in tels], weights=areas)
ref_height = 2180 * u.m
print(f'Reference location: {ref_lon:.6f}°, {ref_lat:.6f}°, {ref_height:}')

print("Transformation using pyproj")
crs = CRS.from_dict(dict(
    proj='tmerc',
    ellps='WGS84',
    lon_0=ref_lon,
    lat_0=ref_lat,
    axis='enu',
    units='m'
))
proj = Transformer.from_crs(crs.geodetic_crs, crs)

for tel in tels:
    proj_easting, proj_northing = proj.transform(tel['loc'].lon.deg, tel['loc'].lat.deg)

    # x_proj is west -> east, y_proj is south -> north.
    # ground frame is x->north, y -> west
    x_proj, y_proj = proj_northing, -proj_easting
    z_proj = (tel['loc'].height - ref_height).to_value(u.m)
    print(f"{tel['name']:10s} x={x_proj: 7.3f} m, y={y_proj: 7.3f} m, z={z_proj: 6.3f} m")

print()
print("Transformation using astropy")
ref = EarthLocation(lon=ref_lon * u.deg, lat=ref_lat * u.deg, height=ref_height)

for tel in tels:
    ground = ground_frame_from_earth_location(tel['loc'], ref)
    print(f"{tel['name']:10s} x={ground.x: 7.3f}, y={ground.y: 7.3f}, z={ground.z: 6.3f}")

Results in (and this is what I would propose for the new locations (the astropy version to be consistent):

Reference location: -17.890879°, 28.761579°, 2180.0 m
Transformation using pyproj
LST-1      x= -5.823 m, y= 60.373 m, z= 19.883 m
MAGIC-1    x= 40.492 m, y=-79.329 m, z= 15.363 m
MAGIC-2    x=-29.833 m, y=-31.180 m, z= 17.364 m

Transformation using astropy
LST-1      x= -5.825 m, y= 60.394 m, z= 19.883 m
MAGIC-1    x= 40.506 m, y=-79.356 m, z= 15.362 m
MAGIC-2    x=-29.843 m, y=-31.191 m, z= 17.364 m

maxnoe commented 1 year ago

The reference location should also be set using the metaparams, which will be used by ctapipe to automatically add this information to the SubarrayDescription.

See e.g. the file sim_telarray/cfg/CTA/CTA-PROD6-metaparam.cfg of the current sim_telarray release:

metaparam global add altitude
#if defined(CTA_NORTH)
   metaparam global set longitude=-17.8920302
   metaparam global set latitude=28.7621661
#elif defined(CTA_SOUTH)
   metaparam global set longitude=-70.316345
   metaparam global set latitude=-24.683429
#endif

maxnoe commented 1 year ago

Ah, and also, relevant xkcd:

[](https://xkcd.com/2170/)

cta-observatory / lst-sim-config

Telescope z positions are raised 10.5 m higher #43