#!/bin/bash

# add: 
# - fkt f dch value? 
# - function for correction factor
# - x-day binning
# < 20121212 data
# zd, th for internal

# todo: check cu-factor for <20121212)

# comparison: 
#  function with select/table: 18.4 sec
#  function with if: 5.7 sec
# function with if: 
#  vorteil: schneller
#  nachteil: keine history
#  nachteil: es braucht 2 fkt (qla, isdc)

# create function CU(night int, ana tinyint) returns double(4,1) deterministic begin declare cu double; set cu=0; select fCU into cu from CU where fValid=1 and fAnalysis=ana and fNight<night order by fNight desc limit 0,1; return cu; end $$

# floor((mjd(timestamp)-mjd(startnight)+0.5)/numdays)
# floor((mjd(timestamp)-mjd(startnight)-0.5)/numdays)


function get_results()
{
   # some basic query parts

   # DataCheck (old)
   ##data with old feedback and/or different bias voltage
   #query=$query" AND fNight>20120420 AND NOT fNight IN (20120406,20120410,20120503) AND"
   ## broken bias channel
   #query=$query" NOT fNight BETWEEN 20121206 AND 20130110"
   # bg-rate cut
   zdparam=" pow(0.753833*cos(Radians(fZenithDistanceMean)), 7.647435)*exp(-5.753686*pow(Radians(fZenithDistanceMean),2.089609))"
   thparam=" pow((if(isnull(fThresholdMinSet),fThresholdMedian,fThresholdMinSet)-329.4203),2)*(-0.0000002044803) "
   param=" (fNumEvtsAfterBgCuts/5-fNumSigEvts)/fOnTimeAfterCuts - "$zdparam" - "$thparam" "
   dchold=" -0.085 < ("$param") "
   dchold=$dchold" AND 0.25 > ("$param") "
   # Datacheck (new) -> combine
   dchval=" fNumEvtsAfterBgCuts/(1.41*POW(fZenithDistanceMean*PI()/180,2)+0.975)/(-7.53e-12*POW(10, LOG10(fThresholdMinSet)*3.69)+1.035)/TIME_TO_SEC(TIMEDIFF(fRunStop,fRunStart))/fEffectiveOn "
   # some semi-automatic datacheck
   dch="  AND (("$dchval" BETWEEN 0.8 AND 1.7 AND fNight BETWEEN 20140520 AND 20150131) " #A
   dch=$dch" OR  ("$dchval" BETWEEN 0.4 AND 1.6 AND fNight BETWEEN 20150201 AND 20150715) " #B
   dch=$dch" OR  ("$dchval" BETWEEN 0.7 AND 1.4 AND fNight BETWEEN 20150716 AND 20160218) " #C
   dch=$dch" OR  ("$dchval" BETWEEN 0.5 AND 1.0 AND fNight > 20160220) " #D
   dch=$dch" OR  ("$dchold" AND fNight<20140520)) " #old

   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
   # some sanity checks
   where=$where" AND fRunTypeKey=1 "
#   where=$where" AND NOT ISNULL(fNumSigEvts) AND NOT ISNULL(fNumBgEvts) "
#   where=$where" AND NOT fRunStart='0000-00-00 00:00:00' AND NOT fRunStop='0000-00-00 00:00:00' "
   # 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
#   crabflux="2.8e-11"
#   crabflux="3.9e-11"
   # range of crab fluxes: 
   # Dortmund: 2.8e-11
   # HESS: 3.37e-11
   # HAWC: 3.01e-11
   # Wue: 3.39e-11 - 3.9e-11 (ISDC analysis)
   # 15-20% difference
   
   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
   
# order information such that it is easily readable by tgraph
# tgraph: X, Y 
# tgraph errors: X, Y, EX, EY

# 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?
# 
# 
# exter	nal
# --------
# 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, 1) 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, 1) 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 "
      #queryint=$queryint" now() as now "
      
      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 "
      #queryext=$queryext" now() as now "
      
      querybase=$from$where
      #if [ $bin -eq -1 ]
      #then
      #   querybase=$querybase" GROUP BY fNight "
      #fi
      #querybase=$querybase" GROUP BY FLOOR((Mjd(fRunStart)-Mjd("$nightmin")-0.5)/"`echo $bin | sed -e 's/-//'`".)"
      querybase=$querybase" GROUP BY num "
      #echo "-----------------"$querybase
      if [ "$ontimelimit" = "" ]
      then 
         querybase=$querybase" HAVING SUM("$ontimeif")>1200 AND NOT SUM(fNumBgEvts)=SUM(fNumSigEvts) ORDER BY num " # 20 min
      else
         querybase=$querybase" HAVING SUM("$ontimeif")>"$ontimelimit" AND NOT SUM(fNumBgEvts)=SUM(fNumSigEvts) 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 "
      #queryint=$queryint" now() as now "
      
      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 "
      #queryext=$queryext" now() as now "
      
      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 NOT sum(o.bgevts)=sum(o.sigevts) AND ontime>0.75*"$bin
      # AND NOT SUM(fNumBgEvts)=SUM(fNumSigEvts)
      
      queryint=$queryint" "$querybase
      queryext=$queryext" "$querybase
   fi
   
   
   fileint=$path"/data/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 root  $dbname -s -e "$queryint" >> $fileint
   #mysql --defaults-file=$sqlpw -u root  $dbname -e "$queryint" 
   
   
   fileext=$path"/data/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 root  $dbname -s -e "$queryext" >> $fileext
   #mysql --defaults-file=$sqlpw -u root  $dbname -e "$queryext"
}

# setup
# db
sqlpw=/home/$USER/.mysql.pw
dbname=factdata20170804
# selection
timeunit=mjd
#bin=20 # min
#bin=0 # period
#bin=-365 # yearly
bin=-1 # nightly
zdmax=90
thmax=1500
path=`dirname $0`
table="AnalysisResultsRunLP"
table="AnalysisResultsAllQLA"
overwrite="yes"
bin=-1 # nightly
nightmin=20111115
nightmax=20171231


# 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



exit

# more examples

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



exit

overwrite="yes"

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

exit

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


exit

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


exit


# 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

exit


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

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


# full sample 
name="Mrk501_all_3d"
source=2
bin=-3
get_results

exit

# full sample 
name="Mrk421_all_3d"
source=1
bin=-3
get_results

name="1959_2016_nightly"
source=7
nightmin=20160301
nightmax=20160831
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

bin=20 # min

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

name="Mrk501_all_20min"
source=2
get_results

name="1959_all_20min"
source=7
get_results

exit

name="Mrk501_all_weekly"
nightmin=20140524
nightmax=20140930
source=2
bin=-7
get_results

exit