source: trunk/DataCheck/Tools/get_data.sh@ 19031

Last change on this file since 19031 was 19031, checked in by Daniela Dorner, 6 years ago
fixed bug
  • Property svn:executable set to *
File size: 14.7 KB
Line 
1#!/bin/bash
2
3# todo
4# - update function for correction
5# - update CU for QLA
6# - add CU for ISDC analysis
7# - add zd, th for internal
8# - add < 20121212 data for QLA
9# - check crab flux
10# - add E2dNdE?
11
12function get_results()
13{
14 # some basic query parts
15
16 # data check based on artificial trigger rate
17 #dch=" AND fR750Cor/fR750Ref >0.93 "
18 dch=" AND fR750Cor/fR750Ref BETWEEN 0.93 AND 1.3 "
19 # ontime
20 ontime1=" TIME_TO_SEC(TIMEDIFF(fRunStop,fRunStart))*fEffectiveOn "
21 ontime2=" fOnTimeAfterCuts "
22 ontimeif=" IF(ISNULL(fEffectiveOn), "$ontime2", "$ontime1") "
23 from=" FROM RunInfo LEFT JOIN "$table" USING (fNight, fRunID) "
24 # time range and source
25 where=" WHERE fSourceKey="$source" AND fNight BETWEEN "$nightmin" AND "$nightmax
26 where=$where" AND NOT ISNULL(fNumExcEvts) "
27 # some sanity checks
28 where=$where" AND fRunTypeKey=1 "
29 # zd cut
30 where=$where" AND fZenithDistanceMax < "$zdmax
31 # th cut
32 where=$where" AND fThresholdMedian < "$thmax
33 where=$where" "$dch
34
35 cufactor=" Avg(25.2) "
36 crabflux="3.37e-11"
37 fluxprec=13
38 crabflux="3.37"
39 fluxprec=2
40
41 case $timeunit in
42 mjd) start=" Mjd(Min(fRunStart)) "
43 stop=" Mjd(MAX(fRunStop)) "
44 deltat=" (Mjd(MAX(fRunStop))-Mjd(Min(fRunStart)))/2 "
45 time=" Mjd(Min(fRunStart))+"$deltat
46 start2=" Mjd(MIN(o.start)) "
47 stop2=" Mjd(MAX(o.stop)) "
48 deltat2=" (Mjd(MAX(o.stop))-Mjd(MIN(o.start)))/2 "
49 time2=" Mjd(MIN(o.start))+"$deltat2
50 ;;
51 unix) start="Unix_timestamp(CONVERT_TZ(Min(fRunStart), '+00:00', 'SYSTEM')) "
52 stop="Unix_timestamp(CONVERT_TZ(Max(fRunStop), '+00:00', 'SYSTEM')) "
53 deltat=" (Unix_timestamp(CONVERT_TZ(Max(fRunStop), '+00:00', 'SYSTEM')) - Unix_timestamp(CONVERT_TZ(Min(fRunStart), '+00:00', 'SYSTEM')))/2 "
54 time=" Unix_timestamp(CONVERT_TZ(Min(fRunStart), '+00:00', 'SYSTEM'))+"$deltat
55 startstop2=" Unix_timestamp(CONVERT_TZ(MIN(o.start), '+00:00', 'SYSTEM')) AS start, "
56 startstop2=$starstop2" Unix_timestamp(CONVERT_TZ(MAX(o.stop), '+00:00', 'SYSTEM')) AS stop, "
57 time2=" (Unix_timestamp(CONVERT_TZ(Max(o.stop), '+00:00', 'SYSTEM')) - Unix_timestamp(CONVERT_TZ(Min(o.start), '+00:00', 'SYSTEM')))/2 "
58 time2=" Unix_timestamp(CONVERT_TZ(Min(o.start), '+00:00', 'SYSTEM'))+"$deltat2
59 ;;
60 *) start=" MIN(fRunStart) "
61 stop=" MAX(fRunStop) "
62 deltat=" sec_to_time(time_to_sec(timediff(MAX(fRunStop), Min(fRunStart)))/2) "
63 time=" addtime(Min(fRunStart), "$deltat") "
64 start2=" MIN(o.start) "
65 stop2=" MAX(o.stop) "
66 deltat2=" sec_to_time(time_to_sec(timediff(MAX(o.stop), Min(o.start)))/2) "
67 time2=" addtime(Min(o.start), "$deltat2") "
68 ;;
69 esac
70 ontime=" SUM("$ontimeif")/60."
71 ontime2=" SUM(o.ot)/60. "
72
73 excrate=" SUM(fNumExcEvts)/SUM("$ontimeif")*3600 "
74 excerr="ExcErr(Sum(fNumSigEvts), SUM(fNumBgEvts))"
75 significance="LiMa(Sum(fNumSigEvts), SUM(fNumBgEvts))"
76 numexc="Sum(fNumExcEvts)"
77 numsig="Sum(fNumSigEvts)"
78 numbg="Sum(fNumBgEvts)"
79 excrateerr=" "$excerr"/SUM("$ontimeif")*3600 "
80 # thomas correction factor
81 correvts=" fNumExcEvts*(pow(cos(fZenithDistanceMean*PI()/180),3)+14.8/21.9*pow(sin(2*fZenithDistanceMean*PI()/180),5))/((1-0.00124/1.21*(if(isnull(fThresholdMinSet),fThresholdMedian,fThresholdMinSet)-500)*(if(isnull(fThresholdMinSet),fThresholdMedian,fThresholdMinSet)>=500))) "
82 correxcrate=" SUM("$correvts")/SUM("$ontimeif")*3600 "
83 # corerr = MMath::ErrorExc(excevtssum+bgevtssum, bgevtssum*5, 0.2)/ontimesum*3600.*corrate/excrate;
84 correxcrateerr=" "$excerr"/SUM("$ontimeif")*3600*SUM("$correvts")/SUM(fNumExcEvts) "
85 # correction on run basis (not relevant for hess)
86 #cu=$correxcrate"/"$cufactor
87 cu=" SUM("$correvts"/CUQLA(fNight))/SUM("$ontimeif")*3600 "
88 #cuerr=$correxcrateerr"/"$cufactor
89 cuerr=" "$excerr"/SUM("$ontimeif")*3600*SUM("$correvts"/CUQLA(fNight))/SUM(fNumExcEvts) "
90 flux=$cu" * "$crabflux
91 fluxerr=$cuerr" * "$crabflux
92
93
94 excrate2=" (SUM(o.sigevts)-SUM(o.bgevts))/SUM(o.ot)*3600 "
95 excerr2="ExcErr(SUM(o.sigevts),SUM(o.bgevts))"
96 significance2="LiMa(SUM(o.sigevts),SUM(o.bgevts))"
97 numexc2="Sum(o.sigevts-o.bgevts)"
98 numsig2="Sum(o.sigevts)"
99 numbg2="Sum(o.bgevts)"
100 excrateerr2=" "$excerr2"/SUM(o.ot)*3600 "
101 correxcrate2=" SUM(o.corevts)/SUM(o.ot)*3600 "
102 correxcrateerr2=" "$excerr2"/SUM(o.ot)*3600*SUM(o.corevts)/(SUM(o.sigevts)-SUM(o.bgevts)) "
103 #cu2=$correxcrate2"/"$cufactor
104 cu2=" SUM(o.corevts/o.cu)/SUM(o.ot)*3600 "
105 #cuerr2=$correxcrateerr2"/"$cufactor
106 cuerr2=" "$excerr2"/SUM(o.ot)*3600*SUM(o.corevts/o.cu)/(SUM(o.sigevts)-SUM(o.bgevts)) "
107 flux2="$cu2*"$crabflux
108 fluxerr2="$cuerr2*"$crabflux
109
110# internal
111# --------
112# timeselect:
113# mjdstar, mjdstop, mjdmean, ontime
114# excselect:
115# excrate, excerr
116# corrected: excrate, excerr
117# CU CUerr
118# flux, fluxerr
119# addselect:
120# signif
121# num exc, num sig, num bg
122# other info: zd? th?
123#
124#
125# external
126# --------
127# time, delta time, start, stop
128# corr-excrate, corr-excerr
129# flux, flux-err
130
131 if [ $bin -le 0 ]
132 then
133 queryint="SELECT "
134 if [ $bin -eq 0 ]
135 then
136 queryint=$queryint" fPeriod as num, "
137 else
138 queryint=$queryint" FLOOR((Mjd(fRunStart)-Mjd("$nightmin")-0.5)/"`echo $bin | sed -e 's/-//'`".) as num, "
139 fi
140 queryint=$queryint" "$time" as time, "$start" as start, "$stop" as stop, "
141 queryint=$queryint" round("$excrate", 1) as excrate, round("$correxcrate", 1) as correxcrate, "
142 queryint=$queryint" round("$cu", 2) as cu, "$flux" as flux, "
143 queryint=$queryint" "$deltat" as deltat, round("$ontime", 1) as ontime, "
144 queryint=$queryint" round("$excrateerr", 1) as excrateerr, round("$correxcrateerr", 1) as correxcrateerr, "
145 queryint=$queryint" round("$cuerr", 2) as cuerr, "$fluxerr" as fluxerr, "
146 queryint=$queryint" round("$significance", 1) as significance, "
147 queryint=$queryint" Min(fNight) as nightmin, "
148 queryint=$queryint" Max(fNight) as nightmax, "
149 queryint=$queryint" "$numexc" as numexc, "
150 queryint=$queryint" "$numsig" as numsig, "
151 queryint=$queryint" "$numbg" as numbg "
152
153 queryext="SELECT "
154 if [ $bin -eq 0 ]
155 then
156 queryext=$queryext" fPeriod as num, "
157 else
158 queryext=$queryext" FLOOR((Mjd(fRunStart)-Mjd("$nightmin")-0.5)/"`echo $bin | sed -e 's/-//'`".) as num, "
159 fi
160 queryext=$queryext" "$time" as time, "$start" as start, "$stop" as stop, "
161 queryext=$queryext" round("$correxcrate", 1) as correxcrate, round("$flux", "$fluxprec") as flux, "
162 queryext=$queryext" "$deltat" as deltat, round("$ontime", 1) as ontime, "
163 queryext=$queryext" round("$correxcrateerr", 1) as correxcrateerr, round("$fluxerr", "$fluxprec") as fluxerr, "
164 queryext=$queryext" round("$significance", 1) as significance "
165
166 querybase=$from$where
167 querybase=$querybase" GROUP BY num "
168 if [ "$ontimelimit" = "" ]
169 then
170 querybase=$querybase" HAVING SUM("$ontimeif")>1200 ORDER BY num " # 20 min
171 else
172 querybase=$querybase" HAVING SUM("$ontimeif")>"$ontimelimit" ORDER BY num "
173 fi
174
175 queryint=$queryint" "$querybase
176 queryext=$queryext" "$querybase
177 else
178 queryint="SELECT "
179 queryint=$queryint" "$time2" as time, "$start2" as start, "$stop2" as stop, "
180 queryint=$queryint" round("$excrate2", 1) as excrate, round("$correxcrate2", 1) as correxcrate, "
181 queryint=$queryint" round("$cu2", 1) as cu, round("$flux2", "$fluxprec") as flux, "
182 queryint=$queryint" round("$excrateerr2", 1) as excrateerr, round("$correxcrateerr2", 1) as correxcrateerr, "
183 queryint=$queryint" "$deltat2" as deltat, round("$ontime2", 1) as ontime, "
184 queryint=$queryint" round("$cuerr2", 1) as cuerr, round("$fluxerr2", "$fluxprec") as fluxerr, "
185 queryint=$queryint" round("$significance2", 1) as significance, "
186 queryint=$queryint" avg(o.night) as night, "
187 queryint=$queryint" "$numexc2" as numexc, "
188 queryint=$queryint" "$numsig2" as numsig, "
189 queryint=$queryint" "$numbg2" as numbg "
190
191 queryext="SELECT "
192 queryext=$queryext" "$time2" as time, "$start2" as start, "$stop2" as stop, "
193 queryext=$queryext" round("$correxcrate2", 1) as correxcrate, round("$flux2", "$fluxprec") as flux, "
194 queryext=$queryext" "$deltat2" as deltat, round("$ontime2", 1) as ontime, "
195 queryext=$queryext" round("$correxcrateerr2", 1) as correxcrateerr, round("$fluxerr2", "$fluxprec") as fluxerr, "
196 queryext=$queryext" round("$significance2", 1) as significance "
197
198 querybase=" FROM (SELECT fNight, @ot:="$ontimeif" AS ot, fRunStart AS start, fRunStop AS stop, fNumSigEvts AS sigevts, fNumBgEvts AS bgevts, "
199 querybase=$querybase" "$correvts" AS corevts, CUQLA(fNight) AS cu, "
200 querybase=$querybase" IF (@night=fNight AND FLOOR((@os+@ot)/"$bin"./60.)<1, @bl, @bl := @bl + 1) AS block, "
201 querybase=$querybase" IF (@night=fNight AND FLOOR((@os+@ot)/"$bin"./60.)<1, @os:=@os + @ot, @os := @ot) AS os, @night :=fNight AS night "
202 querybase=$querybase$from" CROSS JOIN (SELECT @night :=0, @ot :=0, @os :=0, @bl:=0) PARAMS "
203 querybase=$querybase$where" ORDER BY fRunStart) o GROUP BY block HAVING ontime>0.75*"$bin
204
205 queryint=$queryint" "$querybase" order by 'time'"
206 queryext=$queryext" "$querybase" order by 'time'"
207 fi
208
209
210 fileint=$datapath"/FACT_preliminary_"$name"_internal.dat"
211 if [ "$overwrite" = "yes" ]
212 then
213 echo "internal: "$fileint
214 echo "# this file was created at "`date` > $fileint
215 fi
216 if [ $bin -le 0 ]
217 then
218 echo "# numbin time[mjd] start[mjd] stop[mjd] excrate[evts/h] corr.excrate[evts/h] flux[CU] flux[e-11/cm2/s] delta_time[mjd] ontime[min] excrate_err[evts/h] corr.excrate_err[evts/h] flux_err[CU] flux_err[e-11/cm2/s] significance nightmin, nightmax num_exc num_sig num_bg " >> $fileint
219 else
220 echo "# time[mjd] start[mjd] stop[mjd] excrate[evts/h] corr.excrate[evts/h] flux[CU] flux[e-11/cm2/s] delta_time[mjd] ontime[min] excrate_err[evts/h] corr.excrate_err[evts/h] flux_err[CU] flux_err[e-11/cm2/s] significance nightmin, nightmax num_exc num_sig num_bg " >> $fileint
221 fi
222 #echo "$queryint"
223 mysql --defaults-file=$sqlpw -u factread --host=$host $dbname -s -e "$queryint" >> $fileint
224 #mysql --defaults-file=$sqlpw -u factread --host=$host $dbname -e "$queryint"
225
226
227 fileext=$datapath"/FACT_preliminary_"$name".dat"
228 if [ "$overwrite" = "yes" ]
229 then
230 echo "external: "$fileext
231 echo "# this file was created at "`date` > $fileext
232 fi
233 if [ $bin -lt 0 ]
234 then
235 echo "# numbin time[mjd] start[mjd] stop[mjd] corr.excrate[evts/h] flux[e-11/cm2/s] delta_time[mjd] ontime[min] corr.excrate_err[evts/h] flux_err[e-11/cm2/s] significance " >> $fileext
236 else
237 echo "# time[mjd] start[mjd] stop[mjd] corr.excrate[evts/h] flux[e-11/cm2/s] delta_time[mjd] ontime[min] corr.excrate_err[evts/h] flux_err[e-11/cm2/s] significance " >> $fileext
238 fi
239 #echo "$queryext"
240 mysql --defaults-file=$sqlpw -u factread --host=$host $dbname -s -e "$queryext" >> $fileext
241 #mysql --defaults-file=$sqlpw -u factread --host=$host $dbname -e "$queryext"
242}
243
244# setup
245# db
246sqlpw=/home/$USER/.mysql.pw # file with mysql credentials
247#host=lp-fact
248host=10.0.100.21
249#host=localhost
250dbname=factdata # name of database
251# defaults for zd and threshold
252zdmax=90 # all data
253thmax=1500 # all data
254# output path
255path=`dirname $0`
256datapath=$path"/data"
257if ! [ -e $datapath ]
258then
259 mkdir $datapath
260fi
261# time unit
262#timeunit=timestamp # default
263#timeunit=unix
264timeunit=mjd
265# time binning
266# positive values: minutes
267# negative values: days
268# special case 0: period
269# for season binning choose -365 and according start date
270#bin=20 # minutes
271#bin=0 # period
272bin=-1 # nightly
273#bin=-365 # yearly
274# choose analysis
275#table="AnalysisResultsAllQLA" # N/A
276table="AnalysisResultsRunLP" # QLA
277#table="AnalysisResultsRunISDC" # ISDC
278# time range
279nightmin=20111115
280nightmax=20171231
281# overwrite dataset file?
282# (useful to combine different binnings in one file -> set to "no")
283overwrite="yes"
284
285
286# example (adapt to your needs)
287
288# 501 MAGIC
289source=2
290name="Mrk501_2014JulAug"
291bin=-1
292nightmin=20140714
293nightmax=20140805
294get_results
295
296
297
298# end script here
299exit
300
301
302
303#
304# more examples
305#
306
307# Mrk 421
308source=1
309name="Mrk421_nightly"
310bin=-1
311get_results
312name="Mrk421_20min"
313bin=20
314get_results
315name="Mrk421_3d"
316bin=-3
317get_results
318name="Mrk421_10d"
319bin=-10
320get_results
321name="Mrk421_period"
322bin=0
323get_results
324
325
326
327# Mrk 501
328source=2
329name="Mrk501_nightly"
330bin=-1
331get_results
332name="Mrk501_20min"
333bin=20
334get_results
335name="Mrk501_3d"
336bin=-3
337get_results
338name="Mrk501_10d"
339bin=-10
340get_results
341name="Mrk501_period"
342bin=0
343get_results
344
345
346
347# 2344
348source=3
349name="2344_nightly"
350bin=-1
351get_results
352name="2344_20min"
353bin=20
354get_results
355name="2344_period"
356bin=0
357get_results
358
359
360
361# 1959
362source=7
363name="1959_nightly"
364bin=-1
365get_results
366name="1959_20min"
367bin=20
368get_results
369name="1959_period"
370bin=0
371get_results
372
373
374
375# 0323
376source=12
377name="0323_nightly"
378bin=-1
379get_results
380name="0323_20min"
381bin=20
382get_results
383name="0323_period"
384bin=0
385get_results
386
387
388
389# crab
390source=5
391name="Crab_nightly"
392bin=-1
393get_results
394name="Crab_20min"
395bin=20
396get_results
397name="Crab_period"
398bin=0
399get_results
400name="Crab_season"
401bin=-365
402nightmin=20110716
403nightmax=20180716
404get_results
405
406
407
408name="1959_2016"
409source=7
410bin=-1
411nightmin=20160201
412nightmax=20161105
413get_results
414
415name="1959_all_variable"
416overwrite="no"
417source=7
418bin=-365
419nightmin=20120201
420nightmax=20130131
421get_results
422nightmin=20130201
423nightmax=20140131
424get_results
425nightmin=20140201
426nightmax=20150131
427get_results
428bin=0
429nightmin=20150201
430nightmax=20160131
431get_results
432bin=-1
433nightmin=20160201
434nightmax=20170131
435get_results
436bin=0
437nightmin=20170201
438nightmax=20180131
439get_results
440
441
442
443overwrite="yes"
444name="1959_all_variable2"
445overwrite="no"
446source=7
447bin=-365
448nightmin=20120201
449nightmax=20130131
450get_results
451nightmin=20130201
452nightmax=20140131
453get_results
454nightmin=20140201
455nightmax=20150131
456get_results
457bin=0
458nightmin=20150201
459nightmax=20160131
460get_results
461bin=-1
462nightmin=20160201
463nightmax=20160817
464get_results
465bin=0
466nightmin=20160818
467nightmax=20180131
468get_results
469
470
471
472overwrite="yes"
473bin=0
474source=3
475name="2344period"
476get_results
477
478
479
480# flare night (HESS)
481name="Mrk501_10min_flarenight"
482source=2
483bin=10
484nightmin=20140623
485nightmax=20140623
486get_results
487
488
489
490# flare night (HESS)
491name="Mrk501_5min_flarenight"
492source=2
493bin=5
494nightmin=20140623
495nightmax=20140623
496get_results
497
498
499
500
501# full sample
502name="Mrk421_all_nightly"
503source=1
504get_results
505
506name="Mrk501_all_nightly"
507source=2
508get_results
509
510name="1959_all_nightly"
511source=7
512get_results
513
514name="2344_all_nightly"
515source=3
516get_results
517
518
519
520name="HESE20160427"
521source=19
522nightmin=20160425
523bin=-10
524get_results
525
526name="AMON20160731"
527source=21
528nightmin=20160730
529bin=-10
530get_results
531
532
533
Note: See TracBrowser for help on using the repository browser.