fermi-lat / Likelihood

BSD 3-Clause "New" or "Revised" License
2 stars 1 forks source link

Ludicrously high Npred for source out of RoI with energy dispersion #103

Closed jballet closed 3 years ago

jballet commented 3 years ago

I am plagued by a serious bug that completely screws up the results in one of my RoIs. It is a rare bug (only one of my 1500 RoIs is affected) but I don't know how to work around it. Npred is enormous for W51C in the 100-300 MeV band, an extended source 22° off horizontally and 10° off vertically from the center of my RoI, which is a square of 9° half-width. The source is more than 10° off the RoI border, so it is supposed to contribute only a little, but it gets Npred=6 million at edisp_bins=-2, more than everything else in the RoI (including the diffuse emission). I think the problem has something to do with the source map. It is not obviously wrong by itself, but it has two peculiarities:

  1. It stops at half the map (it is zero in the right half), probably because gtsrcmaps considers that it is not worth computing the model so far off
  2. It is non-zero in the first two energy bins (just above 100 MeV) but exactly zero in the next four bins (going higher toward 300 MeV). I suspect this is the root cause of the problem. The extrapolation toward low energies made when edisp_bins is negative probably explodes because it is done in log, or something like this.

The solution may be to change gtsrcmaps to avoid those 0 values, or to change pyLikelihood to teach it how to cope with it. I have simplified the problem to its bare bones and put a test harness at ftp://ftp.cea.fr/incoming/y2k01/jbbglast/SpectralFit_597.tgz (available for one week). The livetime cube is standard and available at /nfs/farm/g/glast/g/catalog/P8_P305. Just launch SpectralFit_597.py. It shows Npred for three values of edisp_bins. The other test sources look all right. I see this both in FT 1.3.5 and 1.4.6.

jballet commented 3 years ago

I have made some progress on this. I now think the bug is in gtsrcmaps. By comparing with another source map for a source closer to the RoI, so that it has no 0 values in any energy bin, and computing what is the offset between the two, I have convinced myself that the first source map (with 237618 in the top left corner) is correct, but the second one (with 36 in the top left corner) is hopelessly wrong. That 36 should be around 120000. So the extrapolation toward low energies of a trend from 36 to 237000 gives rise to very high values, obviously.

donhorner commented 3 years ago

I don't have any help to offer, but I can confirm that the issue also exists in the Fermitools 2.0.8. I ran your test harness, and the value for W51C for edisp_bins=-2 is greater than 6 million. Hopefully, Tom or maybe Eric might have some insight into the underlying issue.

jballet commented 3 years ago

The source map for EXTW51C had been computed some time ago. I have checked this morning that gtsrcmaps in FT 1.4.6 gives the same (probably incorrect) source map, which results in the very large Npred.

eacharles commented 3 years ago

Ok, this is coming from a cut on the padding applied when using a diffuse map as a source. Right now the padding applied is:

double pad_dist =  std::min(10.0,std::max(1.0,meanpsf.containmentRadius(*energy,0.99)));
double radius = std::min(180., data_map_radius + pad_dist);

https://github.com/fermi-lat/Likelihood/blob/50e204a25b6cadaf6b89de364a439fc69acd46fc/src/PSFUtils.cxx#L415

In this particular case, simply changing 0.99 to 0.999 solves the problem, the NPRED for EDISP_BINS =-2, -1, 0 are:

EXTW51C 4.89422300929 EXTW51C 4.32689178795 EXTW51C 3.69384728241

I don't see any reason not to make that change.

jballet commented 3 years ago

I understand that this will make the map size bigger, so it will reduce the part with only 0s. But contrary to what I wrote in my first comment, I am now convinced (from my second comment) that the problem was not the 0s but the wrong (way too low but non-zero) values in the map for the second energy bin. I fear changing the containment will just shift the problem, and it will make the calculation longer for everyone. Did you investigate the origin of those low values in the second energy bin?

eacharles commented 3 years ago

That is the pre convolved map.

On Feb 9, 2021, at 3:30 AM, jballet notifications@github.com wrote:

 I understand that this will make the map size bigger, so it will reduce the part with only 0s. But contrary to what I wrote in my first comment, I am now convinced (from my second comment) that the problem was not the 0s but the wrong (way too low but non-zero) values in the map for the second energy bin. I fear changing the containment will just shift the problem, and it will make the calculation longer for everyone. Did you investigate the origin of those low values in the second energy bin?

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub, or unsubscribe.

eacharles commented 3 years ago

I.e., it will not just change the part with only the zeros but it will soften the edges of the map, decreasing the chances of getting this problem of having a sharply truncated map at one energy bin and not the next.

-e

On Feb 9, 2021, at 3:30 AM, jballet notifications@github.com wrote:

 I understand that this will make the map size bigger, so it will reduce the part with only 0s. But contrary to what I wrote in my first comment, I am now convinced (from my second comment) that the problem was not the 0s but the wrong (way too low but non-zero) values in the map for the second energy bin. I fear changing the containment will just shift the problem, and it will make the calculation longer for everyone. Did you investigate the origin of those low values in the second energy bin?

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub, or unsubscribe.

jballet commented 3 years ago

OK. I just wanted to make sure that you had looked at my second comment.

dagorym commented 3 years ago

You can find a python2.7 version of the FermiTools with this fix as version 1.4.7 the dev channel, i.e.

conda create -n fermi -c conda-forge -c fermi -c fermi/label/dev fermitools=1.4.7

to install.

jballet commented 3 years ago

I have deleted the EXTW51C extension and launched gtsrcmaps from 1.4.7 but I don't see any difference. Could something have gone wrong in the build? Do you see a difference in the source map on your side?

dagorym commented 3 years ago

Yes, it was my mistake not double checking that Eric's change was actually merged into the branch the build used. I'm rebuilding it now. I'll let you know when it is done. It will have the same version number so you'll want to delete the old 1.4.7 and reinstall it. Sorry about that.

On Wed, Feb 10, 2021 at 6:46 AM jballet notifications@github.com wrote:

I have deleted the EXTW51C extension and launched gtsrcmaps from 1.4.7 but I don't see any difference. Could something have gone wrong in the build? Do you see a difference in the source map on your side?

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/fermi-lat/Likelihood/issues/103#issuecomment-776715436, or unsubscribe https://github.com/notifications/unsubscribe-auth/AA3LPSJT5TDQJTPANO4N423S6KEZRANCNFSM4W6RERXQ .

dagorym commented 3 years ago

The rebuild is now done and available.

On Wed, Feb 10, 2021 at 10:13 AM Tom Stephens thomas.stephens@gmail.com wrote:

Yes, it was my mistake not double checking that Eric's change was actually merged into the branch the build used. I'm rebuilding it now. I'll let you know when it is done. It will have the same version number so you'll want to delete the old 1.4.7 and reinstall it. Sorry about that.

On Wed, Feb 10, 2021 at 6:46 AM jballet notifications@github.com wrote:

I have deleted the EXTW51C extension and launched gtsrcmaps from 1.4.7 but I don't see any difference. Could something have gone wrong in the build? Do you see a difference in the source map on your side?

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/fermi-lat/Likelihood/issues/103#issuecomment-776715436, or unsubscribe https://github.com/notifications/unsubscribe-auth/AA3LPSJT5TDQJTPANO4N423S6KEZRANCNFSM4W6RERXQ .

jballet commented 3 years ago

Thanks Tom. The values in the second source map are now reasonable, and the source maps up to 250 MeV contain reasonable and non-zero values as well. Contrary to what I had understood, the sharp drop at about half the RoI is still there, but it is not a problem.

jasercion commented 2 years ago

Can this be closed out?

jballet commented 2 years ago

Hi Joe, Yes, go ahead and close it. Jean

De : jasercion @.> Envoyé : lundi 13 juin 2022 23:12 À : fermi-lat/Likelihood @.> Cc : Ballet Jean @.>; State change @.> Objet : Re: [fermi-lat/Likelihood] Ludicrously high Npred for source out of RoI with energy dispersion (#103)

Can this be closed out?

— Reply to this email directly, view it on GitHubhttps://github.com/fermi-lat/Likelihood/issues/103#issuecomment-1154444077, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AI34LBJWOG7CJHLUPKSZV43VO6P3TANCNFSM4W6RERXQ. You are receiving this because you modified the open/close state.Message ID: @.**@.>>