Skip to content

Conversation

@israelmcmc
Copy link
Contributor

The first and last time bins didn't have a correct exposure due to edge effects and not accounting for the GTI. There were 2 issues:

  1. When using time_ref, the first/last edge could go before/beyond the beginning/end of the TTE time range. This is time for which there is no data but it was still counted as part of the exposure.
  2. The tstart/tstop times corresponded to the first/last event. If you have large statistics this is not really an issue, but once you have low statistics --e.g. you cut in energy first-- then the exposure would be underestimated.

e.g.

from gdt.missions.fermi.gbm.tte import GbmTte
from gdt.core.binning.unbinned import bin_by_time

filename,_ = urllib.request.urlretrieve("https://heasarc.gsfc.nasa.gov/FTP/fermi/data/gbm/triggers/2021/bn210925800/current/glg_tte_n6_bn210925800_v00.fit")
    
tte = GbmTte.open(filename)

phaii = tte.to_phaii(bin_by_time, 0.032, time_ref=0.0)

Before this change would return

$ phaii.data.exposure[:5]
array([0.0319974, 0.031889 , 0.0318864, 0.0318812, 0.0319202])

Where the rate of the first bin is artificially too small

And after:

$ phaii.data.exposure[:5]
array([0.00172335, 0.031889  , 0.0318864 , 0.0318812 , 0.0319202 ])

I know I need to write a test for this and the other open PR. However since these are changes that we'll use for the Burstcube pipeline (we're using the fork https://github.com/BurstCube/gdt-core), I thought it would be better the open the PR sooner than later so people are aware.

@BillCleveland-USRA
Copy link
Contributor

I like how the bounds are enforced with max and min while calculating the exposure, however it should be conditional on the keywords tstart and tstop being included in kwargs.

I'm not sure how the default values for tstart and tstop would be different than the entire time range of the data object if the default arguments are used. Need more discussion on how this routine comes up with the desired answer with the fix.

# Conflicts:
#	src/gdt/core/data_primitives.py
… Adding test to check the default tstart and tstop.
@israelmcmc
Copy link
Contributor Author

Thanks for reviewing, @BillCleveland-USRA. I'm sorry it took me so long to come back to this.

Maybe it's worth decoupling the changes to EventList.bin() (data_primitives.py) and PhotonList.to_phaii() (tte.py)

EventList:
I didn't change the default tstart/tstop here. The only change is that, before, they were not taken into account during the exposure calculation. The code assumed that the data start/stop time matched the first/last bin edges, which is not necessarily true.

PhotonList:
This is where I changed the default. Previously, the implicit default was to use binning method default. In bin_by_time, tstart/tstop are by default the first and last event. This is problematic for the exposure calculation since the start/stop of the observation could have began/stop earlier/later than that first/last photon in the list. For example, if you filter out events based on energy, the time of the first/last event might change, but the exposure shouldn't change.

I'm not sure how the default values for tstart and tstop would be different than the entire time range of the data object if the default arguments are used. Need more discussion on how this routine comes up with the desired answer with the fix.

Yeah, I think I agree. Using the entire time range would be a better default than the GTI. I thought about using the PhotonList.time_range property, but unfortunately it defaults to the time of the first/last event. Should this be changed instead to keep track of the TSTART/TSTOP in the header?

I pushed a couple extra changes:

  • I realized that, in PhotonList.to_phaii(), I should be using obj instead of self. obj is self after cuts.
  • I modified the unit tests for EventList. The reference values also assumed the data fully covered all of the time bins, which is not longer the case if we take into account tstart/tstop.

Some other notes:

  • If we all agree, I can work on changing the behavior of PhotonList.time_range as mention above.
  • I can add unit tests for this change in ``PhotonList.to_phaii` after we resolve the issue above.
  • The following unit tests is failing for me, but it seems unrelated to these changes. Have you seen this before? Any ideas?
tests/core/spectra/test_fitting.py:756 (TestCOBYLA.test_parameters)
self = <tests.core.spectra.test_fitting.TestCOBYLA testMethod=test_parameters>

    def test_parameters(self):
>       self.assertAlmostEqual(self.fitter.parameters[0], 0.05, places=2)
E       AssertionError: np.float64(0.0001000001) != 0.05 within 2 places (np.float64(0.049899999900000006) difference)

# Conflicts:
#	src/gdt/core/tte.py
@israelmcmc
Copy link
Contributor Author

I synced my branch with the main branch and resolved conflicts.

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.

2 participants