ctools

"""Make FITS event lists."""

from glob import glob

import numpy as np

def make_event_list(mssg_max, offset_max):

    """Makes a single FITS event list out of per-run ROOT event lists.

    Takes a few minutes on my laptop to process 20 GB ROOT -> 2 GB FITS.

    TODO: Check if Model++ really used FK5 J2000 coordinates.

    """

    from root_numpy import root2array

    from astropy.table import Table

    from astropy.coordinates import FK5, Galactic

    DIR = '/Users/deil/work/_Data/hess/HESSEVENTLISTS/franzi/WithOwnTables/'

    filenames = glob(DIR + 'EventsList_SpectralHD_Custom_??????.root')[::100]

    selection = '(MeanScaledShowerGoodness_Model < {0})'.format(mssg_max)

    selection += '& (Psi2 < {0})'.format(offset_max ** 2)

    print('Processing {0} files.'.format(len(filenames)))

    array = root2array(filenames, selection=selection)

    print('Selected {0} events.'.format(len(array)))

    print('Making table...')

    table = Table(array)

    # Add header keywords to make the FITS event lists compatible with FTOOLS    

    table.meta['EXTNAME'] = 'EVENTS'

    missing_keywords = ['TSTART', 'TSTOP', 'MJDREF', 'RA_PNT', 'DEC_PNT']

    for missing_keyword in missing_keywords:

        table.meta[missing_keyword] = 0

    print('Computing Galactic coordinates ...')

    ra, dec = table['Lambda'], table['Beta']

    coordinates = FK5(ra, dec, unit=('deg', 'deg'))

    coordinates = coordinates.transform_to(Galactic)

    table['GLON'] = coordinates.lonangle.deg

    table['GLAT'] = coordinates.latangle.deg

    # Add and rename columns to make the FITS event lists compatible with FTOOLS    

    missing_columns = ['EVENT_ID', 'MULTIP', 'DIR_ERR', 'AZ', 'COREX', 'COREY',

                       'CORE_ERR', 'XMAX', 'XMAX_ERR', 'ENERGY_ERR']

    for missing_column in missing_columns:

        table[missing_column] = 0

    table['ALT'] = np.degrees(90 - np.arccos(table['CosZen']))

    col_names = [('MJD', 'TIME'), ('Lambda', 'RA'), ('Beta', 'DEC'),

                 ('NomPosX', 'DETX'), ('NomPosY', 'DETY'), ('Energy', 'ENERGY')]

    for old_name, new_name in col_names:

        table.rename_column(old_name, new_name)

    filename = 'test_events.fits'

    print('Writing file: {0}'.format(filename))

    table.write(filename, overwrite=True)

if __name__ == '__main__':

    make_event_list(mssg_max=-0.6, offset_max=2)