Loads the photon data around a source, subject to a GTI.

The function get_photon_data loads a data dataset of photons in a cone about a specified source, with time, position and energy info, the latter two binned.

These are retrieved from a set of pickled dicts, one per week.

get_photon_data[source]

get_photon_data(config:configuration data, source:Source data, key='')

Parameters:

• source -- PointSource object
• key [''] cache key -- default, use "photons_name", set to None to ignore cache

Steps:

• Read photon data from a Parquet dataset,
• select cone around the source
• use exposure to add exposures
• return DataFrame with columns band time pixel radius

Test reading from a dataset

config = Config(verbose=3)
if config.valid:
    source = PointSource('Geminga')
    print(f'Test loading a photon data set, for source {source.name}')
    photon_data = get_photon_data(config,  source, key=None)
    print(f'Head of the table, length {len(photon_data):,}:\n{photon_data.head()}')
else:
    print('Not testing since no files.')

Test loading a photon data set, for source Geminga
Loading  132 months from Arrow dataset /home/burnett/data/dataset
....................................................................................................................................
	Selected 1,313,726 photons within 5 deg of  (195.13,4.27)
	Energies: 100.0-1000000 MeV
	Dates:    2008-08-04 15:46 - 2019-08-03 01:17
	MJD  :    54682.7          - 58698.1         
Head of the table, length 1,313,726:
   band          time    pixel    radius
0     6  54682.657022  6738278  0.698381
1     3  54682.657934  6761152  2.498099
2     4  54682.658637  6739138  0.290310
3     1  54682.658760  6714890  3.276757
4    11  54682.659099  6736721  4.899003

Parquet conversion code

Copied here, needs integration. Don't know how to make more of those "time_info" filesm, but the data come from the binned FITS files. Need to modify this to read from the FITS files

# #####################################################################################
# #    Parquet code -- TODO just copied, need to test.
# #           from notebooks/code development/parquet_writer.ipynb

# class ParquetConversion(object):
#     import glob, pickle;
#     import healpy
#     import pyarrow as pa
#     import pyarrow.parquet as pq
    
#     def __init__(self, 
#                  data_file_pattern ='$FERMI/data/P8_P305/time_info/month_*.pkl',
#             dataset = '/nfs/farm/g/glast/u/burnett/analysis/lat_timing/data/photon_dataset'):

#         self.files = sorted(glob.glob(os.path.expandvars(data_file_pattern)));
#         print(f'Found {len(self.files)} monthly files with pattern {data_file_pattern}'\
#              f'\nWill store parquet files here: {dataset}')
#         if os.path.exists(dataset):
#             print(f'Dataset folder {dataset} exists')
#         else:
#             os.makedirs(dataset)
#         self.dataset=dataset
            
#     def convert_all(self):
#         files=self.files
#         dataset=self.dataset
#         nside=1024
    
#         def convert(month):

#             infile = files[month-1]
#             print(month, end=',')
#             #print(f'Reading file {os.path.split(infile)[-1]} size {os.path.getsize(infile):,}' )   

#             with open(infile, 'rb') as inp:
#                 t = pickle.load(inp,encoding='latin1')

#             # convert to DataFrame, add month index as new column for partition, then make a Table
#             df = pd.DataFrame(t['timerec'])
#             tstart = t['tstart']
#             df['month']= np.uint8(month)
#             # add a columb with nest indexing -- makes the ring redundant, may remove later
#             df['nest_index'] = healpy.ring2nest(nside, df.hpindex).astype(np.int32)
#             table = pa.Table.from_pandas(df, preserve_index=False)

#             # add to partitioned dataset
#             pq.write_to_dataset(table, root_path=dataset, partition_cols=['month'] )

#         for i in range(len(files)):
#             convert(month=i+1)