Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

h_max wrt array_level #1030

Open
thomasgas opened this issue Mar 25, 2019 · 2 comments
Open

h_max wrt array_level #1030

thomasgas opened this issue Mar 25, 2019 · 2 comments

Comments

@thomasgas
Copy link

The estimation of h_max in HillasReconstructor:

def estimate_h_max(self):

is giving a result wrt the array observation level.
Should we have this number with respect to the sea level?
I used the same piece of code used in the hillas_intersector and inside ImPACT to convert the output of the HillasReconstructor (here called reco_hillas) to the "corrected" one:

def correction_xmax(zenith, height, obs_elev):
    height = height*np.cos(zenith)

    # Add on the height of the detector above sea level
    height = height +  obs_elev
        
    # Lookup this height in the depth tables, the convert Hmax to Xmax
    x_max = thickness_profile(height.to(u.km))
    # Convert to slant depth
    x_max /= np.cos(zenith)
    return x_max

xmax_corrected = correction_xmax(
    np.pi/2 - event.mc.alt.to_value(u.rad), 
    reco_hillas.h_max, 
    event.mcheader.prod_site_alt
)

print(f'reco h_max \t= {reco_hillas.h_max:.2f}')
print(f'mc.energy \t= {event.mc.energy:.2f}')
print(f'from mc \t= {thickness_profile(reco_hillas.h_max.to(u.km)):.2f}')
print(f'event.mc.x_max \t= {event.mc.x_max:.2f}')
print(f'corrected x_max = {xmax_corrected:.2f}')

and this is what i obtain for a huge event, with a good h_max reconstruction (21 telescopes from gamma, Paranal, z20):

reco h_max 	= 9035.94 m
mc.energy 	= 2.60 TeV
from mc 	= 331.01 g / cm2
event.mc.x_max 	= 286.86 g / cm2
corrected x_max = 280.40 g / cm2

We correct it or we leave as it is, specifying that it's wrt the array observation level?

@kosack
Copy link
Contributor

kosack commented Mar 26, 2019

I think sea-level makes the most sense (and is probably what anybody would expect it to be). We should change the Field description to say that it is that explicitly.

This also makes me notice that we need to add the site location into the instrument.SubarrayDescription (not just in the MC header, since it needs to be there for real data files as well)

@thomasgas
Copy link
Author

This goes together with the additional atmospheric profile for La Palma (cta-observatory/ctapipe-extra#24)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants