FT 1.3.1 does not work at all in unbinned mode

[jballet]

The version of unbinnedAnalysis in FT 1.3.1 is completely wrong. It returns minuscule Npred values for point sources, so it finds very large Prefactors, and the fit tends to put most of the counts into the point sources. The attached file contains a simple test harness (testFT_ub.py) illustrating the problem. It writes Npred values for all sources, then fits the parameters (the parameter values can be seen in the Minuit output), then writes Npred and Ts values for all sources. testFT_ub.old is the (normal) output of FT 1.2.23. testFT_ub.new is the (completely wrong) output of FT 1.3.1. testFT_ub.zip

[jballet]

I have noticed that the exposure maps generated by gtexpmaps in FT 1.3.1 are also completely wrong. They look annular instead of filled disks (compare the output below with that in the test harness above). This is not the reason for the problem reported above (the exposure map in the test harness was generated with FT 1.2.23) but it is certainly related. expMap_1001_E19.zip

[dagorym]

Trying to look at this but the test script is referencing a file: $(DATADIR)/Pass8/newdiffuse/InnerGalaxyYB01/isotropic_8YP305_P8R3_SOURCE_V2_v7_YB01.txt that isn't included in the zip file and isn't part of the FermiTools distribution. Where can I find a copy of the file to use for testing?

[jballet]

Here it is, sorry for that. But the problem I see is so enormous that you should see it with whatever isotropic you enter. isotropic_8YP305_P8R3_SOURCE_V2_v7_YB01.txt

[dagorym]

Thanks. You're probably right but its easier to be able to exactly reproduce the situation you had so as not to get distracted by other variables.

[dagorym]

Interestingly, this works just fine in the current 1.9.9 version that uses Python 3 and ROOT 6. I get the exact same results as the old 1.2.23 version of the science tools. So that means it's either something very recent that isn't in that branch but is in the 1.3.1 version, or is some strange interaction between the dependencies that had to be shoe-horned to get ROOT 6 to work with Python 2.7.

I'm leaning toward the latter simply because, IIRC, I pulled all the changes made to create the 1.3.1 version into the Python 3 update branch shortly after the merge of the ROOT 6 change into master that went into the 1.3.1 version of the tools.

I'll keep looking to try to track down the exact place in the history where the change happened.

[dagorym]

This bug was introduced by the fix to this issue: https://github.com/fermi-lat/Likelihood/issues/70 that fixed the gtexpmap failure. Changing to code back to the old check (ie. psi<roi_radius instead of cos(roi_radius-psi)<0.99) makes this issue go away but reintroduces the previous bug.

I think the real question here is the tolerance. The old gtexpmap bug was introduced when the psi value was allowed to go right up to the edge of the roi_radius. the limit of 0.99 means that it's not allowed to get to within 8.1 degrees. If we move the tolerance to 0.999 that allows the psi value to get to within 2.56 degrees of the roi radius and at 0.9999 we get to 0.81 degrees.

At 0.999 this test still fails and the gtexpmap test passes. At 0.9999 both of these tests past.

So, from a software perspective, with just these two tests, there is a way to have them both pass. Whether or not one value over the other makes sense scientifically, is a different question. For now. I'll commit the change setting the limit to 0.9999. If we find another test case where this fails we can revisit it.

[jballet]

What do you mean with "this test passes"? That you get exactly the same Npred and Ts values for all sources, and the same exposure map, in the development branch as in the official FT 1.2.23? If so, this looks good to me.

[dagorym]

I didn't try to create an exposure map to compare but the numerical output of the test script you supplied is exactly the same now.

[jballet]

I must admit I don't understand how this integration works, but 0.81 deg from the ROI boundary is not very small. No source in the test harness is within that distance of the RoI border, and the isotropic probably builds on the exposure map, which was computed beforehand. So I am not sure it is testing the interesting region. I have prepared a different test harness in which I have manually placed two sources very close to the border (504F-0256 inside and S8008-0794 outside). If that gives the same output as in 1.2.23, all will be OK. testFT_border.zip I am not sure gtexpmap uses the same part of the code, so compare the exposure map (built by calling gtexpmap.sh) as well, to see whether anything changes within 0.8° of the border.

[dagorym]

Where is the gtexpmap.sh file? It wasn't in the testFT_border.zip file. Was it supposed to be?

[jballet]

I thought I had put it in the expMap_1001_E19.zip file. Here it is (renamed it .txt because .sh is not accepted) gtexpmap.txt

[dagorym]

Thanks.

[dagorym]

Jean, exactly which values should be matching in these runs? There are a lot of moving parts here and I'm try to assess what should be stable. I'm guessing that it's the NPred and TS printed at the very end for each source that is the most important although the former doesn't change before and after the fit so maybe it's not as relevant?

The change introduced by Eric to fix the other issue has the effect, I believe, of slightly shrinking the size of the region evaluated in the ROI by limiting how close to the edge you can get before things are ignored. Because of this, I think things are going to be slightly different no matter what happens if the sources are close to the edge of the ROI.

[dagorym]

Okay, regardless of the answers to the above questions. I think I've figured this out.

As I decreased the region of exclusion (basically adding more 9's, i.e. .9999 -> .99999999) things got closer to the original value but then started bouncing around a bit. In fact, using mum <1 which should have been equivalent to the previous psi< roi_radius, gave a slightly different value for the S8008 source, while all the others stayed the same, which prompted the previous questions. (The previous test that Eric was trying to fix also failed.)

However, as I was writing that previous comment, I realized that if psi is slightly greater than roi_radius, the cosine of the difference would be the same as if it was slightly smaller. This means that events outside the radius could be slipping in resulting in those different values. So I kept in the mum check to guard against the previous error that the change was introduced to fix and also added back in the psi<roi_radius check (via a logical AND). That fixed the discrepancy and the new version of the code exactly matches the older version of the code for both test cases here and the gtexpmap test that Eric was fixing also passes.

I'll look at the expmap generation script that Jean sent and verify that those are the same now as well.

[jballet]

The numbers to check are the Npred for the two sources that I noted above (504F-0256 and S8008-0794). They are fixed because there are no events there so of course they cannot be fitted. I could even have removed the fit from the script and called Npred directly. TS is irrelevant. Adding back psi < roi_radius is definitely a good idea. Our PSF gets down to 0.03° (PSF3 at high energy), so this is the minimum width of the ramp at the RoI border. Rather than checking the cosinus (which is closer to 1 than 1E-7 at a distance of 0.01°) it may be simpler to replace psi < roi_radius by psi < roi_radius - 1E-5, say (assuming psi and roi_radius are expressed in radians).

[dagorym]

Okay. comparing an exposure map made with the fixed version of the tools and the exposure map Jean supplied as part of the test case (made with the 1.2.23 version of the tools) shows no differences (other than a different order of header keywords). The data is the same in both.

I think this issue is now resolved.

[jballet]

Excellent! All clear then.