source: branches/trigger_burst_research/calculate_burst_ratio_file.py@ 18679

Last change on this file since 18679 was 18304, checked in by dneise, 9 years ago
gave output file, more useful name
File size: 4.3 KB
Line 
1# coding: utf-8
2import progressbar
3import calendar
4import numpy as np
5from astropy.io import fits
6import pandas as pd
7from sqlalchemy import create_engine
8
9night_int = 20150721
10
11
12def fjd(datetime_inst):
13 """ convert datetime instance to FJD
14 """
15 return calendar.timegm(datetime_inst.utctimetuple()) / (24.*3600.)
16
17def create_DB_connection():
18 from ConfigParser import SafeConfigParser
19 config = SafeConfigParser()
20 config.optionxform = str # this make the parsing case sensitive
21 config.read('config.ini')
22
23 factdb = create_engine(
24 "mysql+mysqldb://{user}:{pw}@{host}/{db}".format(
25 user=config.get('database', 'user'),
26 pw=config.get('database', 'password'),
27 host=config.get('database', 'host'),
28 db=config.get('database', 'database'),
29 )
30 )
31
32 return factdb
33
34
35def get_all_data_runs_from_run_db(factdb):
36 keys = [
37 'fRunID',
38 'fNight',
39 'fRunStart',
40 'fRunStop',
41 'fZenithDistanceMean',
42 'fAzimuthMean',
43 ]
44
45 sql_query = """SELECT {comma_sep_keys}
46 FROM RunInfo WHERE fRunTypeKey=1 ORDER BY fNight;
47 """
48 sql_query = sql_query.format(
49 comma_sep_keys=', '.join(keys),
50 )
51
52 data_runs = pd.read_sql_query(
53 sql_query,
54 factdb,
55 parse_dates=['fRunStart', 'fRunStop'],
56 )
57
58 return data_runs
59
60
61
62def get_trigger_rates(night_int, base_path='/fact/aux/'):
63 night_string = str(night_int)
64 fits_file = fits.open(
65 base_path+'{y}/{m}/{d}/{n}.FTM_CONTROL_TRIGGER_RATES.fits'.format(
66 n=night_string,
67 y=night_string[0:4],
68 m=night_string[4:6],
69 d=night_string[6:8],
70 )
71 )
72 trigger_rates = fits_file[1].data
73 return trigger_rates
74
75
76def add_ratio_and_more_to_dataframe(data_runs):
77
78 # create new columns and pre-assign some hopefully good NULL-like-values
79 data_runs['number_of_measurements'] = 0
80 data_runs['median_of_board_rates'] = np.nan
81 data_runs['std_dev_of_board_rates'] = np.nan
82 data_runs['fBoardTriggerRateRatioAboveThreshold'] = np.nan
83
84 trigger_rates = None
85 last_night = None
86
87 progress = progressbar.ProgressBar(widgets=[progressbar.Bar('=', '[', ']'), ' ',
88 progressbar.Percentage(), ' ',
89 progressbar.ETA()])
90
91 for dataframe_index in progress(data_runs.index):
92 this_run = data_runs.ix[dataframe_index]
93 if last_night != this_run['fNight']:
94 try:
95 trigger_rates = get_trigger_rates(this_run['fNight'])
96 except (IOError, ValueError):
97 continue
98 last_night = this_run['fNight']
99
100 try:
101 mask = (
102 (trigger_rates['Time'] > fjd(this_run['fRunStart'])) *
103 (trigger_rates['Time'] < fjd(this_run['fRunStop']))
104 )
105 except ValueError:
106 continue
107
108 if mask.sum() == 0:
109 continue
110 try:
111 this_board_rates = trigger_rates['BoardRate'][mask]
112 except KeyError:
113 continue
114
115 N = this_board_rates.shape[0]
116 data_runs.ix[dataframe_index, 'number_of_measurements'] = N
117
118 left, med, right = np.percentile(
119 this_board_rates,
120 [50-20, 50, 50+20])
121 something_like_width = (right - left)/2.
122 # 40% of the area of the normal distrubution
123 # fall between -0.52*sigma and +0.52*sigma.
124 # The estimation of sigma in a distorted CDF
125 # is less affected near the median. (but of course less accurate as well.)
126 std_dev = something_like_width / 0.52
127
128 median_board_rates = np.median(this_board_rates, axis=1)
129 over = (median_board_rates > med+3*std_dev).sum()
130
131 data_runs.ix[dataframe_index, 'median_of_board_rates'] = med
132 data_runs.ix[dataframe_index, 'std_dev_of_board_rates'] = std_dev
133 data_runs.ix[dataframe_index, 'fBoardTriggerRateRatioAboveThreshold'] = float(over)/N
134
135def main(night_int):
136 factdb = create_DB_connection()
137 data_runs = get_all_data_runs_from_run_db(factdb)
138 add_ratio_and_more_to_dataframe(data_runs)
139 return data_runs
140
141if __name__ == '__main__':
142 data_runs = main(night_int)
143 data_runs.to_csv('burst_ratio.csv')
Note: See TracBrowser for help on using the repository browser.