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

weighting MC gamma to MAGIC log-parabola spectrum #302

Merged
merged 12 commits into from
Mar 30, 2020

Conversation

SeiyaNozaki
Copy link
Collaborator

Here I'd like to implement the code for weighting MC gamma to LogParabola spectrum.
Main difference is:

  • can select spectral shape for weighting, "PowerLaw" or "LogParabola" in mc/mc.py("rate", "weight")
  • apply LogParabola in mc/sensitivity.py("find_best_cuts_sensitivity", "sensitivity")
    Any comments are welcome.

lstchain/mc/mc.py Show resolved Hide resolved

return rate


def weight(emin, emax, sim_sp_idx, rate, nev, w_param):
def weight(shape, emin, emax, sim_sp_idx, rate, nev, w_param):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above, explain w_param, especially signs.

Comment on lines 357 to 359
elif(shape == "LogParabola"):
dFdE, crab_par = crab_magic(energy)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

General comment: for differential sensitivity calculations, using the HEGRA power-law spectrum, or a more accurate log-parabola parametrization, will make no difference. The effect would just go through the difference in hardness within each energy bin , and hence is very small for the usually narrow (5 bins/dex) bins we use.

The log-parabola weighting is instead crucial for comparisons of Crab data with MC gammas, for example to obtain parameter distributions, because in that case the energy range spanned by events entering one such comparison may be arbitrarily wide.

@SeiyaNozaki
Copy link
Collaborator Author

Thanks for the comments @moralejo !

I see, then here I'd like to implement the change for only 'rate', 'weight' function, not sensitivity codes. (If needed in the future, I'll pull request again).
I'll modify the codes following your comments (add explanation of parameters).

@rlopezcoto
Copy link
Contributor

Hi @SeiyaNozaki, I agree with @moralejo that the log-parabola will not make a big difference for the differential sensitivity calculations, but this will be useful for other calculations, thanks for implementing it!
Could you maybe think on adding some test unit along with this PR?

@SeiyaNozaki
Copy link
Collaborator Author

SeiyaNozaki commented Mar 17, 2020

Hi @rlopezcoto Sorry I forgot to answer your message.

Could you maybe think on adding some test unit along with this PR?

For now, no. Just I wanted to implement this function for anyone to use log-parabola weighting.

param should include 'f0','e0','alpha'
dFdE = f0 * np.power(E / e0, alpha)

if shape is 'LogParabola':
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @SeiyaNozaki, I'm happy with the changes and would propose this to be accepted if @moralejo approves, just one additional question: if the shape given is LogParabola and no beta parameter is given, will the call throw an error or is there any default value?
Also, what if shape == PowerLaw and I introduce a beta parameter in the call?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @rlopezcoto ,
For the first case(LogParabola without beta), we will get error.
But the second case(PowerLaw with beta), no error and just other parameters are used for the calcurlation.
Would it be better to add the function to check if the param match the shape?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what kind of error in the first case? (a generic python one I guess).
It would be great if you could check that we are giving the correct parameters depending on the function introduced.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what kind of error in the first case? (a generic python one I guess).

Yes, just a KeyError.

>>> crab_par_hegra
{'f0': <Quantity 2.83e-11 1 / (cm2 s TeV)>, 'alpha': -2.62, 'e0': <Quantity 1. TeV>}
>>> rate_ = rate("LogParabola", 0.1 * u.TeV, 100 * u.TeV, crab_par_hegra, 0, 0)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/cta/Work/Install/anaconda3/envs/cta-dev/lib/python3.7/site-packages/lstchain-0.4.5.post9+git66dd663-py3.7.egg/lstchain/mc/mc.py", line 107, in rate
    log_parabola =  LogParabolaSpectralModel.from_log10(amplitude=param['f0'], reference=param['e0'], alpha=-1*param['alpha'], beta=-1*param['beta'])
KeyError: 'beta'

It would be great if you could check that we are giving the correct parameters depending on the function introduced.

I'll introduce such a function in a few days.

@SeiyaNozaki
Copy link
Collaborator Author

I added the function to check if the input parameters correspond to the expected function.

Copy link
Contributor

@rlopezcoto rlopezcoto left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks a lot @SeiyaNozaki!
If @moralejo approves, I'll merge

@SeiyaNozaki
Copy link
Collaborator Author

Sorry, I got an unit error with Log Parabola function, let me check a little...

@SeiyaNozaki
Copy link
Collaborator Author

I fixed the unit error, sorry.
Now I'm completely sure that this function works fine!

@rlopezcoto rlopezcoto merged commit d83b098 into cta-observatory:master Mar 30, 2020
@SeiyaNozaki SeiyaNozaki deleted the weight_to_magic branch April 29, 2020 16:30
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants