Overview
Fermi-LAT weekly data files are extracted from the GSFC FTP server ,
with subfolders for the photon data, photon and spacecraft data, spacecraft. It is described here
The class FermiData downloads these to temporary files and constructs a dict for each week with
contents
photons: a table, one entry per selected photon with columns, converted with
get_ft1_datarun_ref(uint8), index into therunslist- time since the run start, in 2 $\mu$s intervals (uint32)
- energy and event type (uint8)
- position as HEALPix index with nside, nest from Config, currently 1024, True (uint32)
sc_data: a table, one entry per 30-s interval, with columns, all float32, converted with
get_ft2_info- start/stop time
- S/C direction
- zenith direction
- runs list of run numbers
- gti_times: an array of interleaved start/stop intervals
- file_date: modification date for the FT1 file at GSFC.
These dict objects are saved in a folder specfied by Config().datapath. Its contents can be
checked with check_data, and updated using update_data. The contents of a week can be visually
inspected with plot_week.
Run numbers and timing
The run number is an integer correspnding to the MET at the start of the run. For a week, with 15 orbits/day, we expect ~105 runs. We save a table of the run numbers per week, and use a uint8 index in the photon table.
A run is typically an orbit or less, at most 6 ks. Integerizing the offset from the run start, or the run number, using 32 bits, one has 5e5 intervals/s, so we choose 2$\mu$s. Thus we generate the time for each event from three sources: the 2$\mu$s time offset for the run, the index of the run number into the runs table, and the runs table itself.
GTI
GTI, for Good Time Interval, is a mechanism to select event times, and corresponding spacecraft information if exposure is needed. There is a section for this information in the FT1 photon data files, but, for the GSFC-generated weekly files we use here, we do not use it, since the photon data itself respects it, as do the space craft information in the FT2 files. This code is intended to support external use.
display_markdown("""\
## FT1 and FT2 -- data access
Fermi data is distributed in two FITS file formats:
- FT1: photon data
- FT2: spacecraft history
""")
display_markdown( """### Test download from GSFC and data handling
Download, then process FT1 and FT2 files.
""")
gd = GSFCweekly(Config(verbose=3))
if gd.config.valid:
ft2, ft1 = gd.download(-2)
v1 = get_ft1_data(gd.config, Path(gd.local_path)/ft1);
print(f'Contents of FT1 summary output dict\n'
' key value type')
for (k,x) in v1.items():
print(f' {k:10} {x.__class__.__name__}')
v2 = get_ft2_info(gd.config, Path(gd.local_path)/ft2)
print(f'Contents of FT2 summary output dict\n'
' key value type')
for (k,x) in v2.items():
print(f' {k:10} {x.__class__.__name__}')
gti = v1['gti_times']
gti-gti[0]
g = np.diff(gti)[1::2]
b = np.diff(gti)[0::2]
gti[-1]-gti[0],pd.Series(g).describe(), pd.Series(b).describe(), np.sum(g), np.sum(b)
display_markdown("""### Test GTI
Extract photon data for a week containing two solar flares, check that
they get excluded.
""")
interval = (55990, 55995)
dv = DataView(interval, gti=GTI('nosolarflares.gti'))
print(dv)
wk =dv.photon_week(dv[0])
t = wk['time']
dv_nocut = DataView(interval)
wk_nocut = dv_nocut.photon_week(dv_nocut[0])
bins = np.linspace(t[0],t[-1], 100)
fig, ax = plt.subplots(figsize=(6,2))
ax.hist(wk_nocut['time'], bins=bins, color='lightgrey',label='all data');
ax.hist(t,bins=bins, color='green', label='selected');
ax.set(xlabel='MJD', xlim=interval)
ax.legend();
debug
display_markdown("""## PixelView
Manage the data as a FITS file accepted by pointlike
This class uses a DataView, or a FITS file to load a set of pixels with counts
for a time interval.
The FITS format has two HDUs with columns as shown.
__BANDS__
* NSIDE (J) -- power of 2 up to 4096
* E_MIN (D), E_MAX (D)
* EVENT_TYPE (J) -- 0 or 1 for Front, Back
__SKYMAP__
* PIX (I) -- pixel index (RING) ordering now, but can be specified in header.
* CHANNEL (I) -- band index
* VALUE (J) -- number of events
Currently the wtlike photons are all stored with nside 1024, while the pointlike files
go up to 4096, simply because that it the largest for 32 bits.
""")
pmap = pm.get_map(6,64)
pm.mollview(8, 32,max=4)
#pm.mollview(6); # second time works--mystery TODO
bands.dtype
nside = bands.NSIDE
e_min = bands.E_MIN
e_max = bands.E_MAX
dr3pfile = Path(os.path.expandvars(
'$FERMI/pointlike-data/12years_zmax100_4bpd_v2.fits'))
assert dr3pfile.is_file()
hdus = fits.open(dr3pfile)
hdus.info()