source: branches/trigger_burst_research/calculate_burst_ratio_file.py@ 18296

# coding: utf-8
import progressbar
import calendar
import numpy as np
from astropy.io import fits
import pandas as pd
from sqlalchemy import create_engine

night_int = 20150721


def fjd(datetime_inst):
    """ convert datetime instance to FJD
    """
    return calendar.timegm(datetime_inst.utctimetuple()) / (24.*3600.)

def create_DB_connection():
    from ConfigParser import SafeConfigParser
    config = SafeConfigParser()
    config.optionxform = str  # this make the parsing case sensitive
    config.read('config.ini')

    factdb = create_engine(
        "mysql+mysqldb://{user}:{pw}@{host}/{db}".format(
            user=config.get('database', 'user'),
            pw=config.get('database', 'password'),
            host=config.get('database', 'host'),
            db=config.get('database', 'database'),
        )
    )

    return factdb


def get_all_data_runs_from_run_db(factdb):
    keys = [
        'fRunID',
        'fNight',
        'fRunStart',
        'fRunStop',
        'fZenithDistanceMean',
        'fAzimuthMean',
        ]

    sql_query = """SELECT {comma_sep_keys}
        FROM RunInfo WHERE fRunTypeKey=1 ORDER BY fNight;
        """
    sql_query = sql_query.format(
        comma_sep_keys=', '.join(keys),
        )
    
    data_runs = pd.read_sql_query(
        sql_query,
        factdb,
        parse_dates=['fRunStart', 'fRunStop'],
    )

    return data_runs



def get_trigger_rates(night_int, base_path='/fact/aux/'):
    night_string = str(night_int)
    fits_file = fits.open(
        base_path+'{y}/{m}/{d}/{n}.FTM_CONTROL_TRIGGER_RATES.fits'.format(
            n=night_string,
            y=night_string[0:4],
            m=night_string[4:6],
            d=night_string[6:8],
            )
        )
    trigger_rates = fits_file[1].data
    return trigger_rates


def add_ratio_and_more_to_dataframe(data_runs):

    # create new columns and pre-assign some hopefully good NULL-like-values
    data_runs['number_of_measurements'] = 0
    data_runs['median_of_board_rates'] = np.nan
    data_runs['std_dev_of_board_rates'] = np.nan
    data_runs['fBoardTriggerRateRatioAboveThreshold'] = np.nan

    trigger_rates = None
    last_night = None

    progress = progressbar.ProgressBar(widgets=[progressbar.Bar('=', '[', ']'), ' ',
                                            progressbar.Percentage(), ' ',
                                            progressbar.ETA()])

    for dataframe_index in progress(data_runs.index):
        this_run = data_runs.ix[dataframe_index]
        if last_night != this_run['fNight']:
            try:
                trigger_rates = get_trigger_rates(this_run['fNight'])
            except (IOError, ValueError):
                continue
            last_night = this_run['fNight']

        try:
            mask = (
                (trigger_rates['Time'] > fjd(this_run['fRunStart'])) *
                (trigger_rates['Time'] < fjd(this_run['fRunStop']))
                )
        except ValueError:
            continue
        
        if mask.sum() == 0:
            continue
        try:
            this_board_rates = trigger_rates['BoardRate'][mask]
        except KeyError:
            continue

        N = this_board_rates.shape[0]
        data_runs.ix[dataframe_index, 'number_of_measurements'] = N
        
        left, med, right = np.percentile(
            this_board_rates,
            [50-20, 50, 50+20])
        something_like_width = (right - left)/2.
        # 40% of the area of the normal distrubution
        # fall between -0.52*sigma and +0.52*sigma.
        # The estimation of sigma in a distorted CDF 
        # is less affected near the median. (but of course less accurate as well.)
        std_dev = something_like_width / 0.52
        
        median_board_rates = np.median(this_board_rates, axis=1)
        over = (median_board_rates > med+3*std_dev).sum()

        data_runs.ix[dataframe_index, 'median_of_board_rates'] = med
        data_runs.ix[dataframe_index, 'std_dev_of_board_rates'] = std_dev
        data_runs.ix[dataframe_index, 'fBoardTriggerRateRatioAboveThreshold'] = float(over)/N

def main(night_int):
    factdb = create_DB_connection()
    data_runs = get_all_data_runs_from_run_db(factdb)
    add_ratio_and_more_to_dataframe(data_runs)
    return data_runs

if __name__ == '__main__':
    data_runs = main(night_int)
    data_runs.to_csv('foo.bar')