#!/bin/bash

# todo
# - update function for correction
# - update CU for QLA
# - add CU for ISDC analysis
# - add zd, th for internal
# - add < 20121212 data for QLA
# - check crab flux
# - add E2dNdE?

function get_results()
{
   # some basic query parts

   # data check based on artificial trigger rate
   #dch=" AND fR750Cor/fR750Ref >0.93 "
   dch=" AND fR750Cor/fR750Ref BETWEEN 0.93 AND 1.3 "
   # ontime 
   ontime1=" TIME_TO_SEC(TIMEDIFF(fRunStop,fRunStart))*fEffectiveOn "
   ontime2=" fOnTimeAfterCuts "
   ontimeif=" IF(ISNULL(fEffectiveOn), "$ontime2", "$ontime1") "
   from=" FROM RunInfo LEFT JOIN "$table" USING (fNight, fRunID) "
   # time range and source
   where=" WHERE fSourceKey="$source" AND fNight BETWEEN "$nightmin" AND "$nightmax
   where=$where" AND NOT ISNULL(fNumExcEvts) "
   # some sanity checks
   where=$where" AND fRunTypeKey=1 "
   # zd cut
   where=$where" AND fZenithDistanceMax < "$zdmax
   # th cut
   where=$where" AND fThresholdMedian < "$thmax
   where=$where" "$dch
   
   cufactor=" Avg(25.2) "
   crabflux="3.37e-11"
   fluxprec=13
   crabflux="3.37"
   fluxprec=2
   
   case $timeunit in
      mjd)  start=" Mjd(Min(fRunStart)) "
            stop=" Mjd(MAX(fRunStop)) "
            deltat=" (Mjd(MAX(fRunStop))-Mjd(Min(fRunStart)))/2 "
            time=" Mjd(Min(fRunStart))+"$deltat
            start2=" Mjd(MIN(o.start)) "
            stop2=" Mjd(MAX(o.stop)) "
            deltat2=" (Mjd(MAX(o.stop))-Mjd(MIN(o.start)))/2 "
            time2=" Mjd(MIN(o.start))+"$deltat2
            ;;
      unix) start="Unix_timestamp(CONVERT_TZ(Min(fRunStart), '+00:00', 'SYSTEM')) "
            stop="Unix_timestamp(CONVERT_TZ(Max(fRunStop), '+00:00', 'SYSTEM')) "
            deltat=" (Unix_timestamp(CONVERT_TZ(Max(fRunStop), '+00:00', 'SYSTEM')) - Unix_timestamp(CONVERT_TZ(Min(fRunStart), '+00:00', 'SYSTEM')))/2 "
            time="  Unix_timestamp(CONVERT_TZ(Min(fRunStart), '+00:00', 'SYSTEM'))+"$deltat
            startstop2=" Unix_timestamp(CONVERT_TZ(MIN(o.start), '+00:00', 'SYSTEM')) AS start, "
            startstop2=$starstop2" Unix_timestamp(CONVERT_TZ(MAX(o.stop), '+00:00', 'SYSTEM')) AS stop, "
            time2="  (Unix_timestamp(CONVERT_TZ(Max(o.stop), '+00:00', 'SYSTEM')) - Unix_timestamp(CONVERT_TZ(Min(o.start), '+00:00', 'SYSTEM')))/2 "
            time2="  Unix_timestamp(CONVERT_TZ(Min(o.start), '+00:00', 'SYSTEM'))+"$deltat2
            ;;
      *)    start=" MIN(fRunStart) "
            stop=" MAX(fRunStop) "
            deltat=" sec_to_time(time_to_sec(timediff(MAX(fRunStop), Min(fRunStart)))/2) "
            time=" addtime(Min(fRunStart), "$deltat") "
            start2=" MIN(o.start) "
            stop2=" MAX(o.stop) "
            deltat2=" sec_to_time(time_to_sec(timediff(MAX(o.stop), Min(o.start)))/2) "
            time2=" addtime(Min(o.start), "$deltat2") "
            ;;
   esac
   ontime=" SUM("$ontimeif")/60."
   ontime2=" SUM(o.ot)/60. "
   
   excrate=" SUM(fNumExcEvts)/SUM("$ontimeif")*3600 "
   excerr="ExcErr(Sum(fNumSigEvts), SUM(fNumBgEvts))"
   significance="LiMa(Sum(fNumSigEvts), SUM(fNumBgEvts))"
   numexc="Sum(fNumExcEvts)"
   numsig="Sum(fNumSigEvts)"
   numbg="Sum(fNumBgEvts)"
   excrateerr=" "$excerr"/SUM("$ontimeif")*3600 "
   # thomas correction factor
   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))) "
   correxcrate=" SUM("$correvts")/SUM("$ontimeif")*3600 "
   # corerr = MMath::ErrorExc(excevtssum+bgevtssum, bgevtssum*5, 0.2)/ontimesum*3600.*corrate/excrate;
   correxcrateerr=" "$excerr"/SUM("$ontimeif")*3600*SUM("$correvts")/SUM(fNumExcEvts) "
   # correction on run basis (not relevant for hess)
   #cu=$correxcrate"/"$cufactor
   cu=" SUM("$correvts"/CUQLA(fNight))/SUM("$ontimeif")*3600 "
   #cuerr=$correxcrateerr"/"$cufactor
   cuerr=" "$excerr"/SUM("$ontimeif")*3600*SUM("$correvts"/CUQLA(fNight))/SUM(fNumExcEvts) "
   flux=$cu" * "$crabflux
   fluxerr=$cuerr" * "$crabflux

   
   excrate2=" (SUM(o.sigevts)-SUM(o.bgevts))/SUM(o.ot)*3600 "
   excerr2="ExcErr(SUM(o.sigevts),SUM(o.bgevts))"
   significance2="LiMa(SUM(o.sigevts),SUM(o.bgevts))"
   numexc2="Sum(o.sigevts-o.bgevts)"
   numsig2="Sum(o.sigevts)"
   numbg2="Sum(o.bgevts)"
   excrateerr2=" "$excerr2"/SUM(o.ot)*3600 "
   correxcrate2=" SUM(o.corevts)/SUM(o.ot)*3600 "
   correxcrateerr2=" "$excerr2"/SUM(o.ot)*3600*SUM(o.corevts)/(SUM(o.sigevts)-SUM(o.bgevts)) "
   #cu2=$correxcrate2"/"$cufactor
   cu2=" SUM(o.corevts/o.cu)/SUM(o.ot)*3600 "
   #cuerr2=$correxcrateerr2"/"$cufactor
   cuerr2=" "$excerr2"/SUM(o.ot)*3600*SUM(o.corevts/o.cu)/(SUM(o.sigevts)-SUM(o.bgevts)) "
   flux2="$cu2*"$crabflux
   fluxerr2="$cuerr2*"$crabflux

# internal 
# --------
# timeselect:
#  mjdstar, mjdstop, mjdmean, ontime
# excselect:
#  excrate, excerr
#  corrected: excrate, excerr
#  CU CUerr
#  flux, fluxerr
# addselect: 
#  signif
#  num exc, num sig, num bg
# other info: zd? th?
# 
# 
# external
# --------
# time, delta time, start, stop
# corr-excrate, corr-excerr
# flux, flux-err
   
   if [ $bin -le 0 ]
   then 
      queryint="SELECT "
      if [ $bin -eq 0 ]
      then 
         queryint=$queryint" fPeriod as num, "
      else
         queryint=$queryint" FLOOR((Mjd(fRunStart)-Mjd("$nightmin")-0.5)/"`echo $bin | sed -e 's/-//'`".) as num, "
      fi
      queryint=$queryint" "$time" as time, "$start" as start, "$stop" as stop, "
      queryint=$queryint" round("$excrate", 1) as excrate, round("$correxcrate", 1) as correxcrate, "
      queryint=$queryint" round("$cu", 2) as cu, "$flux" as flux, "
      queryint=$queryint" "$deltat" as deltat, round("$ontime", 1) as ontime, "
      queryint=$queryint" round("$excrateerr", 1) as excrateerr, round("$correxcrateerr", 1) as correxcrateerr, "
      queryint=$queryint" round("$cuerr", 2) as cuerr, "$fluxerr" as fluxerr, "
      queryint=$queryint" round("$significance", 1) as significance, "
      queryint=$queryint" Min(fNight) as nightmin, "
      queryint=$queryint" Max(fNight) as nightmax, "
      queryint=$queryint" "$numexc" as numexc, "
      queryint=$queryint" "$numsig" as numsig, "
      queryint=$queryint" "$numbg" as numbg "
      
      queryext="SELECT "
      if [ $bin -eq 0 ]
      then 
         queryext=$queryext" fPeriod as num, "
      else
         queryext=$queryext" FLOOR((Mjd(fRunStart)-Mjd("$nightmin")-0.5)/"`echo $bin | sed -e 's/-//'`".) as num, "
      fi
      queryext=$queryext" "$time" as time, "$start" as start, "$stop" as stop, "
      queryext=$queryext" round("$correxcrate", 1) as correxcrate, round("$flux", "$fluxprec") as flux, "
      queryext=$queryext" "$deltat" as deltat, round("$ontime", 1) as ontime, "
      queryext=$queryext" round("$correxcrateerr", 1) as correxcrateerr, round("$fluxerr", "$fluxprec") as fluxerr, "
      queryext=$queryext" round("$significance", 1) as significance "
      
      querybase=$from$where
      querybase=$querybase" GROUP BY num "
      if [ "$ontimelimit" = "" ]
      then 
         querybase=$querybase" HAVING SUM("$ontimeif")>1200 ORDER BY num " # 20 min
      else
         querybase=$querybase" HAVING SUM("$ontimeif")>"$ontimelimit" ORDER BY num "
      fi
      
      queryint=$queryint" "$querybase
      queryext=$queryext" "$querybase
   else
      queryint="SELECT "
      queryint=$queryint" "$time2" as time, "$start2" as start, "$stop2" as stop, "
      queryint=$queryint" round("$excrate2", 1) as excrate, round("$correxcrate2", 1) as correxcrate, "
      queryint=$queryint" round("$cu2", 1) as cu, round("$flux2", "$fluxprec") as flux, "
      queryint=$queryint" round("$excrateerr2", 1) as excrateerr, round("$correxcrateerr2", 1) as correxcrateerr, "
      queryint=$queryint" "$deltat2" as deltat, round("$ontime2", 1) as ontime, "
      queryint=$queryint" round("$cuerr2", 1) as cuerr, round("$fluxerr2", "$fluxprec") as fluxerr, "
      queryint=$queryint" round("$significance2", 1) as significance, "
      queryint=$queryint" avg(o.night) as night, "
      queryint=$queryint" "$numexc2" as numexc, "
      queryint=$queryint" "$numsig2" as numsig, "
      queryint=$queryint" "$numbg2" as numbg "
      
      queryext="SELECT "
      queryext=$queryext" "$time2" as time, "$start2" as start, "$stop2" as stop, "
      queryext=$queryext" round("$correxcrate2", 1) as correxcrate, round("$flux2", "$fluxprec") as flux, "
      queryext=$queryext" "$deltat2" as deltat, round("$ontime2", 1) as ontime, "
      queryext=$queryext" round("$correxcrateerr2", 1) as correxcrateerr, round("$fluxerr2", "$fluxprec") as fluxerr, "
      queryext=$queryext" round("$significance2", 1) as significance "
      
      querybase=" FROM (SELECT fNight, @ot:="$ontimeif" AS ot, fRunStart AS start, fRunStop AS stop, fNumSigEvts AS sigevts, fNumBgEvts AS bgevts, "
      querybase=$querybase" "$correvts" AS corevts, CUQLA(fNight) AS cu, "
      querybase=$querybase" IF (@night=fNight AND FLOOR((@os+@ot)/"$bin"./60.)<1, @bl, @bl := @bl + 1) AS block, "
      querybase=$querybase" IF (@night=fNight AND FLOOR((@os+@ot)/"$bin"./60.)<1, @os:=@os + @ot, @os := @ot) AS os, @night :=fNight AS night "
      querybase=$querybase$from" CROSS JOIN (SELECT @night :=0, @ot :=0, @os :=0, @bl:=0) PARAMS "
      querybase=$querybase$where" ORDER BY fRunStart) o GROUP BY block HAVING ontime>0.75*"$bin
      
      queryint=$queryint" "$querybase" order by 'time'"
      queryext=$queryext" "$querybase" order by 'time'"
   fi
   
   
   fileint=$datapath"/FACT_preliminary_"$name"_internal.dat"
   if [ "$overwrite" = "yes" ]
   then 
      echo "internal: "$fileint
      echo "# this file was created at "`date` > $fileint
   fi
   if [ $bin -le 0 ]
   then 
      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
   else
      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
   fi
   #echo "$queryint"
   mysql --defaults-file=$sqlpw -u factread --host=$host $dbname -s -e "$queryint" >> $fileint
   #mysql --defaults-file=$sqlpw -u factread --host=$host $dbname -e "$queryint" 
   
   
   fileext=$datapath"/FACT_preliminary_"$name".dat"
   if [ "$overwrite" = "yes" ]
   then 
      echo "external: "$fileext
      echo "# this file was created at "`date` > $fileext
   fi
   if [ $bin -lt 0 ]
   then 
      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
   else
      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
   fi
   #echo "$queryext"
   mysql --defaults-file=$sqlpw -u factread --host=$host $dbname -s -e "$queryext" >> $fileext
   #mysql --defaults-file=$sqlpw -u factread --host=$host $dbname -e "$queryext"
}

# setup
# db
sqlpw=/home/$USER/.mysql.pw # file with mysql credentials
#host=lp-fact
host=10.0.100.21
#host=localhost
dbname=factdata # name of database
# defaults for zd and threshold
zdmax=90 # all data
thmax=1500 # all data
# output path
path=`dirname $0`
datapath=$path"/data"
if ! [ -e $datapath ]
then
   mkdir $datapath
fi
# time unit 
#timeunit=timestamp # default
#timeunit=unix
timeunit=mjd
# time binning 
# positive values: minutes
# negative values: days
# special case 0: period
# for season binning choose -365 and according start date
#bin=20 # minutes
#bin=0 # period
bin=-1 # nightly
#bin=-365 # yearly
# choose analysis 
#table="AnalysisResultsAllQLA" # N/A
table="AnalysisResultsRunLP" # QLA 
#table="AnalysisResultsRunISDC"  # ISDC
# time range
nightmin=20111115
nightmax=20171231
# overwrite dataset file? 
# (useful to combine different binnings in one file -> set to "no")
overwrite="yes"


# example (adapt to your needs)

# 501 MAGIC 
source=2
name="Mrk501_2014JulAug"
bin=-1
nightmin=20140714
nightmax=20140805
get_results



# end script here
exit



#
# more examples
#

# Mrk 421
source=1
name="Mrk421_nightly"
bin=-1
get_results
name="Mrk421_20min"
bin=20
get_results
name="Mrk421_3d"
bin=-3
get_results
name="Mrk421_10d"
bin=-10
get_results
name="Mrk421_period"
bin=0
get_results



# Mrk 501
source=2
name="Mrk501_nightly"
bin=-1
get_results
name="Mrk501_20min"
bin=20
get_results
name="Mrk501_3d"
bin=-3
get_results
name="Mrk501_10d"
bin=-10
get_results
name="Mrk501_period"
bin=0
get_results



# 2344
source=3
name="2344_nightly"
bin=-1
get_results
name="2344_20min"
bin=20
get_results
name="2344_period"
bin=0
get_results



# 1959
source=7
name="1959_nightly"
bin=-1
get_results
name="1959_20min"
bin=20
get_results
name="1959_period"
bin=0
get_results



# 0323
source=12
name="0323_nightly"
bin=-1
get_results
name="0323_20min"
bin=20
get_results
name="0323_period"
bin=0
get_results



# crab
source=5
name="Crab_nightly"
bin=-1
get_results
name="Crab_20min"
bin=20
get_results
name="Crab_period"
bin=0
get_results
name="Crab_season"
bin=-365
nightmin=20110716
nightmax=20180716
get_results



name="1959_2016"
source=7
bin=-1
nightmin=20160201
nightmax=20161105
get_results

name="1959_all_variable"
overwrite="no"
source=7
bin=-365
nightmin=20120201
nightmax=20130131
get_results
nightmin=20130201
nightmax=20140131
get_results
nightmin=20140201
nightmax=20150131
get_results
bin=0
nightmin=20150201
nightmax=20160131
get_results
bin=-1
nightmin=20160201
nightmax=20170131
get_results
bin=0
nightmin=20170201
nightmax=20180131
get_results



overwrite="yes"
name="1959_all_variable2"
overwrite="no"
source=7
bin=-365
nightmin=20120201
nightmax=20130131
get_results
nightmin=20130201
nightmax=20140131
get_results
nightmin=20140201
nightmax=20150131
get_results
bin=0
nightmin=20150201
nightmax=20160131
get_results
bin=-1
nightmin=20160201
nightmax=20160817
get_results
bin=0	
nightmin=20160818
nightmax=20180131
get_results



overwrite="yes"
bin=0
source=3
name="2344period"
get_results



# flare night (HESS)
name="Mrk501_10min_flarenight"
source=2
bin=10
nightmin=20140623
nightmax=20140623
get_results



# flare night (HESS)
name="Mrk501_5min_flarenight"
source=2
bin=5
nightmin=20140623
nightmax=20140623
get_results




# full sample 
name="Mrk421_all_nightly"
source=1
get_results

name="Mrk501_all_nightly"
source=2
get_results

name="1959_all_nightly"
source=7
get_results

name="2344_all_nightly"
source=3
get_results



name="HESE20160427"
source=19
nightmin=20160425
bin=-10
get_results

name="AMON20160731"
source=21
nightmin=20160730
bin=-10
get_results