Ignore:
Timestamp:
10/07/19 12:57:31 (5 years ago)
Author:
Daniela Dorner
Message:
fixed errors for corrected LCs, added new options (cut on correction-factors, cut for rows with empty values), fixed typo in header
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/DataCheck/Tools/get_data.sh

    r19570 r19731  
    137137         echo "#   Calima Cut: dust < "$dust" ug/cm3"
    138138      fi
     139      if [ "$factorcut" != "" ]
     140      then
     141         echo "#   Factor Cut: zdfactor*thfactor > "$factorcut
     142      fi
    139143      if [ "$usedch" == "yes" ]
    140144      then
     
    175179      zdfactor="(pow(cos(fZenithDistanceMean*PI()/180),3)+14.8/21.9*pow(sin(2*fZenithDistanceMean*PI()/180),5))"
    176180      #zdfactor="(1/(pow(cos(fZenithDistanceMean*PI()/180),3)+14.8/21.9*pow(sin(2*fZenithDistanceMean*PI()/180),5)))"
    177       thfactor="(1-0.00124/1.21*(fThresholdMinSet-500)*(fThresholdMinSet>=500))"                                 
     181      thfactor="(1-0.00124/1.21*(fThresholdMinSet-500)*(fThresholdMinSet>=500))"
    178182      # ETh: 1100 GeV
    179183      crabflux="1.81"
     
    207211   cufactor="fCU"
    208212   # some calculations
    209    excerr="ExcErr(Sum(fNumSigEvts), SUM(fNumBgEvts))"
     213   #excerr="ExcErr(SUM(fNumSigEvts), SUM(fNumBgEvts))"
     214   excerr="ExcErr(fNumSigEvts, fNumBgEvts)"
    210215   CU="SUM(fNumExcEvts/"$cufactor")/SUM("$ontimeif")*3600"
    211    CUerr=$excerr"/SUM("$ontimeif")*3600*SUM(fNumExcEvts/"$cufactor")/SUM(fNumExcEvts)"
     216   #CUerr=$excerr"/SUM("$ontimeif")*3600*SUM(fNumExcEvts/"$cufactor")/SUM(fNumExcEvts)"
     217   CUerr="SQRT(SUM(POW("$excerr"/"$cufactor",2)))/SUM("$ontimeif")*3600"
    212218   CUcor="SUM("$correvts"/"$cufactor")/SUM("$ontimeif")*3600"
    213    CUcorerr=$excerr"/SUM("$ontimeif")*3600*SUM("$correvts"/"$cufactor")/SUM(fNumExcEvts)"
    214    excerr2="ExcErr(SUM(o.sigevts),SUM(o.bgevts))"
     219   #CUcorerr=$excerr"/SUM("$ontimeif")*3600*SUM("$correvts"/"$cufactor")/SUM(fNumExcEvts)"
     220   CUcorerr="SQRT(SUM(POW("$excerr"/"$cufactor"/"$zdfactor"/"$thfactor",2)))/SUM("$ontimeif")*3600"
     221   #excerr2="ExcErr(SUM(o.sigevts),SUM(o.bgevts))"
     222   excerr2="ExcErr(o.sigevts,o.bgevts)"
    215223   CU2="SUM((o.sigevts-o.bgevts)/o.cufactor)/SUM(o.ot)*3600"
    216    CUerr2=$excerr2"/SUM(o.ot)*3600*SUM((o.sigevts-o.bgevts)/o.cufactor)/(SUM(o.sigevts)-SUM(o.bgevts))"
    217    CUcor2="SUM(o.corevts/o.cufactor)/SUM(o.ot)*3600"
    218    CUcorerr2=$excerr2"/SUM(o.ot)*3600*SUM(o.corevts/o.cufactor)/(SUM(o.sigevts)-SUM(o.bgevts))"
     224   #CUerr2=$excerr2"/SUM(o.ot)*3600*SUM((o.sigevts-o.bgevts)/o.cufactor)/(SUM(o.sigevts)-SUM(o.bgevts))"
     225   CUerr2="SQRT(SUM(POW("$excerr2"/o.cufactor,2)))/SUM(o.ot)*3600"
     226   #CUcor2="SUM(o.corevts/o.cufactor)/SUM(o.ot)*3600"
     227   CUcor2="SUM(o.excevts/o.zdfactor/o.thfactor/o.cufactor)/SUM(o.ot)*3600"
     228   #CUcorerr2=$excerr2"/SUM(o.ot)*3600*SUM(o.corevts/o.cufactor)/(SUM(o.sigevts)-SUM(o.bgevts))"
     229   CUcorerr2="SQRT(SUM(POW("$excerr2"/o.cufactor/o.zdfactor/o.thfactor,2)))/SUM(o.ot)*3600"
     230   
     231   #"ROUND(SQRT(SUM(POW(ExcErr(fNumSigEvts, fNumBgEvts)/pow(cos(fZenithDistanceMean*PI()/180)*exp(1-cos(fZenithDistanceMean*PI()/180)),4.5)/(1.37-IF(ISNULL(fThresholdMinSet),fThresholdMedian,fThresholdMinSet)*0.00118), 2)))/SUM( IF(ISNULL(fEffectiveOn),  fOnTimeAfterCuts ,  TIME_TO_SEC(TIMEDIFF(fRunStop,fRunStart))*fEffectiveOn ) )*3600, 1)"
    219232   
    220233   # columns to be selected
     
    223236   excrate=" ROUND(SUM(fNumExcEvts)/SUM("$ontimeif")*3600, 1) AS excrate"
    224237   significance="ROUND(LiMa(Sum(fNumSigEvts), SUM(fNumBgEvts)), 1) AS significance"
    225    numexc="Sum(fNumExcEvts) AS numexc"
    226    numsig="Sum(fNumSigEvts) AS numsig"
    227    numbg="Sum(fNumBgEvts) AS numbg"
    228    excrateerr=" ROUND("$excerr"/SUM("$ontimeif")*3600, 1) AS excrateerr"
     238   numexc="SUM(fNumExcEvts) AS numexc"
     239   numsig="SUM(fNumSigEvts) AS numsig"
     240   numbg="SUM(fNumBgEvts) AS numbg"
     241   #excrateerr=" ROUND("$excerr"/SUM("$ontimeif")*3600, 1) AS excrateerr"
     242   excrateerr=" ROUND(ExcErr(SUM(fNumSigEvts), SUM(fNumBgEvts))/SUM("$ontimeif")*3600, 1) AS excrateerr"
    229243   correxcrate=" ROUND(SUM("$correvts")/SUM("$ontimeif")*3600, 1) AS correxcrate"
    230    correxcrateerr=" ROUND("$excerr"/SUM("$ontimeif")*3600*SUM("$correvts")/SUM(fNumExcEvts), 1) AS correxcrateerr"
     244   #correxcrateerr=" ROUND("$excerr"/SUM("$ontimeif")*3600*SUM("$correvts")/SUM(fNumExcEvts), 1) AS correxcrateerr"
     245   correxcrateerr=" ROUND(SQRT(SUM(POW("$excerr"/"$zdfactor"/"$thfactor", 2)))/SUM("$ontimeif")*3600, 1) AS correxcrateerr"
    231246   cu=" ROUND("$CU", 2) AS cu"
    232247   cuerr=" ROUND("$CUerr", 2) AS cuerr"
    233248   cucor=" ROUND("$CUcor", 2) AS cucor"
    234    cucorerr=" ROUND("$CUcorerr", 2) AS cucorerr"
     249   cucorerr=" ROUND("$CUcorerr", 2) AS 'cucorerr'"
    235250   flux="ROUND("$CU" * "$crabflux", 2) AS flux"
    236251   fluxerr="ROUND("$CUerr" * "$crabflux", 2) AS fluxerr"
     
    239254   # for minute binning
    240255   ontime2=" ROUND(SUM(o.ot)/60., 1) AS ontime"
    241    excrate2=" ROUND((SUM(o.sigevts)-SUM(o.bgevts))/SUM(o.ot)*3600, 1) AS excrate"
     256   #excrate2=" ROUND((SUM(o.sigevts)-SUM(o.bgevts))/SUM(o.ot)*3600, 1) AS excrate"
     257   excrate2=" ROUND((SUM(o.excevts))/SUM(o.ot)*3600, 1) AS excrate"
    242258   significance2=" ROUND(LiMa(SUM(o.sigevts),SUM(o.bgevts)), 1) AS significance"
    243    numexc2="Sum(o.sigevts-o.bgevts) AS numexc"
    244    numsig2="Sum(o.sigevts) AS numsig"
    245    numbg2="Sum(o.bgevts) AS numbg"
    246    excrateerr2=" ROUND("$excerr2"/SUM(o.ot)*3600, 1) AS excrateerr"
    247    correxcrate2=" ROUND(SUM(o.corevts)/SUM(o.ot)*3600, 1) AS correxcrate"
    248    correxcrateerr2=" ROUND("$excerr2"/SUM(o.ot)*3600*SUM(o.corevts)/(SUM(o.sigevts)-SUM(o.bgevts)), 1) AS correxcrateerr"
     259   #numexc2="Sum(o.sigevts-o.bgevts) AS numexc"
     260   numexc2="SUM(o.excevts) AS numexc"
     261   numsig2="SUM(o.sigevts) AS numsig"
     262   numbg2="SUM(o.bgevts) AS numbg"
     263   excrateerr2=" ROUND(ExcErr(SUM(o.sigevts),SUM(o.bgevts))/SUM(o.ot)*3600, 1) AS excrateerr"
     264   correxcrate2=" ROUND(SUM(o.excevts/o.zdfactor/o.thfactor)/SUM(o.ot)*3600, 1) AS correxcrate"
     265   #correxcrateerr2=" ROUND("$excerr2"/SUM(o.ot)*3600*SUM(o.corevts)/(SUM(o.sigevts)-SUM(o.bgevts)), 1) AS correxcrateerr"
     266   correxcrateerr2=" ROUND(SQRT(SUM(POW("$excerr2"/o.zdfactor/o.thfactor,2)))/SUM(o.ot)*3600,1) AS correxcrateerr"
    249267   cu2=" ROUND("$CU2", 2) AS cu"
    250268   cuerr2=" ROUND("$CUerr2", 2) AS cuerr"
    251269   cucor2=" ROUND("$CUcor2", 2) AS cucor"
    252    cucorerr2=" ROUND("$CUcorerr2", 2) AS cucorerr"
     270   cucorerr2=" ROUND("$CUcorerr2", 2) AS 'cucorerr'"
    253271   flux2="ROUND("$CU2" * "$crabflux", "$fluxprec") AS flux"
    254272   fluxerr2="ROUND("$CUerr2" *"$crabflux", "$fluxprec") AS fluxerr"
     
    306324   # some sanity checks
    307325   where=$where" AND fRunTypeKey=1 "
     326   # remove empty rows
     327   if [ "$rmemptyrows" == "yes" ]
     328   then
     329      where=$where" AND NOT ISNULL(fZenithDistanceMean) AND (NOT ISNULL(fThresholdMinSet) OR NOT ISNULL(fThresholdMedian)) AND NOT ISNULL(fCU) "
     330   fi
    308331   # zd cut
    309332   if [ "$zdmax" != "" ]
     
    334357      where=$where" AND "$lightcut
    335358   fi
     359#   if [ "$factorcut" != "" ]
     360#   then
     361#      where=$where" AND ( ("$zdfactor") * ("$thfactor")) > "$factorcut" "
     362#   fi
    336363   querybase=$from$where
    337364
     
    358385      if [ $bin -eq 0 ]
    359386      then
    360          orderby=" fPeriod "
     387         orderby=" ORDER BY fPeriod "
    361388         #querystart=$querystart" fPeriod AS num, "
    362389         queryend=" GROUP BY fPeriod "
    363390      else
    364391         num=" FLOOR((Mjd(fRunStart)-Mjd("$nightmin")-0.5)/"`echo $bin | sed -e 's/-//'`".) "
    365          orderby=$num
     392         orderby=" ORDER BY "$num
    366393         #querystart=$querystart" FLOOR((Mjd(fRunStart)-Mjd("$nightmin")-0.5)/"`echo $bin | sed -e 's/-//'`".) AS num, "
    367394         queryend=" GROUP BY "$num
     
    370397      if [ "$ontimelimit" = "" ]
    371398      then
    372          queryend=$queryend" HAVING SUM("$ontimeif")>1200 ORDER BY "$orderby
     399         queryend=$queryend" HAVING SUM("$ontimeif")>1200 "
    373400      else
    374          queryend=$queryend" HAVING SUM("$ontimeif")>"$ontimelimit" ORDER BY "$orderby
     401         queryend=$queryend" HAVING SUM("$ontimeif")>"$ontimelimit
    375402      fi
     403#      if [ "$ontimelimit" = "" ]
     404#      then
     405#         queryend=$queryend" HAVING abs(cucorerr) < 1 AND SUM("$ontimeif")>1200 ORDER BY "$orderby
     406#      else
     407#         queryend=$queryend" HAVING abs(cucorerr) < 1 AND SUM("$ontimeif")>"$ontimelimit" ORDER BY "$orderby
     408#      fi
    376409     
    377410      # internal
     
    385418      queryint=$queryint" MIN("$zenith"Min) AS zdmin, MAX("$zenith"Max) AS zdmax, "
    386419      queryint=$queryint" MIN("$thresh") AS thmin, MAX("$thresh") AS thmax, "
    387       queryint=$queryint" ROUND(AVG("$cufactor"), 1) AS cufactor, ROUND(AVG(fR750Cor), 2) AS R750cor,  ROUND(AVG(fR750Ref), 2) AS R750ref "
    388       queryint=$queryint" "$querybase" "$querydch" "$queryend
     420      queryint=$queryint" ROUND(AVG("$cufactor"), 1) AS cufactor, ROUND(AVG("$zdfactor"), 1) AS zdfactor, ROUND(AVG("$thfactor"), 1) AS thfactor, "
     421      queryint=$queryint" ROUND(AVG(fR750Cor), 2) AS R750cor,  ROUND(AVG(fR750Ref), 2) AS R750ref "
     422      if [ "$factorcut" == "" ]
     423      then
     424         queryint=$queryint" "$querybase" "$querydch" "$queryend" "$orderby
     425      else
     426#         queryint=$queryint" "$querybase" "$querydch" "$queryend" "$orderby
     427         queryint=$queryint" "$querybase" "$querydch" "$queryend" AND ( zdfactor * thfactor > "$factorcut") "$orderby
     428      fi
    389429     
    390430      # for collaborators
     
    394434      querycol=$querycol" "$excrateerr", "$correxcrateerr", "$cucorerr", "$fluxcorerr", "
    395435      querycol=$querycol" "$significance
    396       querycol=$querycol" "$querybase" "$querydch" "$queryend
     436      if [ "$factorcut" == "" ]
     437      then
     438         querycol=$querycol" "$querybase" "$querydch" "$queryend" "$orderby
     439      else
     440#         querycol=$querycol" "$querybase" "$querydch" "$queryend" "$orderby
     441         querycol=$querycol" "$querybase" "$querydch" "$queryend" AND (AVG("$zdfactor") * AVG("$thfactor") > "$factorcut") "$orderby
     442      fi
    397443     
    398444      # external
    399445      # no datacheck applied for external files
    400       queryext=$querystart" "$excrate", "$deltat", "$excrateerr" "$querybase" "$queryend
     446      # no cut on factors
     447      queryext=$querystart" "$excrate", "$deltat", "$excrateerr" "$querybase" "$queryend" "$orderby
    401448     
    402449   else
     
    408455      querybase=" FROM (SELECT fNight, fZenithDistanceMin AS zdmin, fZenithDistanceMax AS zdmax, "$thresh" AS th, "
    409456      querybase=$querybase" fR750Cor AS R750cor, fR750Ref AS R750ref, "$cufactor" AS cufactor, "
     457      querybase=$querybase" "$zdfactor" AS zdfactor, "$thfactor" AS thfactor, "
    410458      querybase=$querybase" @ot:="$ontimeif" AS ot, fRunStart AS start, fRunStop AS stop, "
    411       querybase=$querybase" fNumSigEvts AS sigevts, fNumBgEvts AS bgevts, "$correvts" AS corevts, "
     459      #querybase=$querybase" fNumSigEvts AS sigevts, fNumBgEvts AS bgevts, "$correvts" AS corevts, "
     460      querybase=$querybase" fNumSigEvts AS sigevts, fNumBgEvts AS bgevts, fNumExcEvts AS excevts, "
    412461      querybase=$querybase" IF (@night=fNight AND FLOOR((@os+@ot)/"$bin"./60.)<1, @bl, @bl := @bl + 1) AS block, "
    413462      querybase=$querybase" IF (@night=fNight AND FLOOR((@os+@ot)/"$bin"./60.)<1, @os:=@os + @ot, @os := @ot) AS os, @night :=fNight AS night "
     
    421470      queryint=$queryint" "$excrateerr2", "$cuerr2", "$fluxerr2", "$correxcrateerr2", "$cucorerr2", "$fluxcorerr2", "
    422471      queryint=$queryint" "$significance2", "
    423       queryint=$queryint" avg(o.night) AS night, "
     472      #queryint=$queryint" avg(o.night) AS night, "
     473      queryint=$queryint" min(o.night) AS nightmin, max(o.night) AS nightmax, "
    424474      queryint=$queryint" "$numexc2", "$numsig2", "$numbg2", "
    425475      queryint=$queryint" MIN(o.zdmin) AS zdmin, MAX(o.zdmax) AS zdmax, MIN(o.th) AS thmin, MAX(o.th) AS thmax, "
     
    480530   echo "#" >> $fileint
    481531   headerint="# time["$timeunit"] start["$timeunit"] stop["$timeunit"] excess-rate[evts/h] flux[CU] flux[e-11ph/cm2/s] corrected_excess-rate[evts/h] corrected_flux[CU] corrected_flux[e-11ph/cm2/s] (stop-start)/2["$timeunit"] ontime[min]"
    482    headerint=$headerint" excess-rate_error[evts/h] flux_error[CU] flux_error[e-11ph/cm2/s] corrected_excess-rate_error[evts/h] corrected_flux_error[CU] corrected_flux_error[e-11/cm2/s] significance night num_exc num_sig num_bg "
    483    headerint=$headerint" zdmin zdmax thmin thmax avg(cufactor) avg(R750cor) avg(R750ref) "
     532   headerint=$headerint" excess-rate_error[evts/h] flux_error[CU] flux_error[e-11ph/cm2/s] corrected_excess-rate_error[evts/h] corrected_flux_error[CU] corrected_flux_error[e-11/cm2/s] significance nightmin nightmax num_exc num_sig num_bg "
     533   headerint=$headerint" zdmin zdmax thmin thmax avg(cufactor) avg(zdfactor) avg(thfactor) avg(R750cor) avg(R750ref) "
    484534   echo $headerint >> $fileint
    485535   #echo "$queryint"
     
    639689# use different conversion from CU to fluxes
    640690#crabfluxconv="2.5"
     691# remove lines which contain NULL (e.g. zd-corrected flux when zd=NULL)
     692rmemptyrows="yes"
     693# remove lines with zdfactor * thfactor > value
     694factorcut=0.3
    641695
    642696
     
    655709# -------------------------------------------------------------------------------------- #
    656710
    657 # some test
    658 table="AnalysisResultsRunCutsLC" # CutsLC
    659 bin=0
     711
     712datapath=/home/dorner/FACT.analysis/corrected.lcs.for.collaborators/lcs
     713
     714# put your data request here, examples below
     715
     716# LCs for 421 campaign with Astrosat Jan 2019
     717table="AnalysisResultsRunLP"
     718# Mrk 421
     719source=1
     720# nightly
     721bin=-1
     722nightmin=20190109
     723nightmax=20190113
     724name="QLA_Mrk421_nightly_Astrosat-Jan-2019"
     725get_results
     726# 20min
     727bin=20
     728name="QLA_Mrk421_20min_Astrosat-Jan-2019"
     729get_results
     730# 30min
     731bin=30
     732name="QLA_Mrk421_30min_Astrosat-Jan-2019"
     733get_results
     734# 40min
     735bin=40
     736name="QLA_Mrk421_40min_Astrosat-Jan-2019"
     737get_results
     738# 60min
     739bin=60
     740name="QLA_Mrk421_60min_Astrosat-Jan-2019"
     741get_results
     742# 80min
     743bin=80
     744name="QLA_Mrk421_80min_Astrosat-Jan-2019"
     745get_results
     746# 90min
     747bin=90
     748name="QLA_Mrk421_90min_Astrosat-Jan-2019"
     749get_results
     750
     751
     752
     753
     754exit
     755
     756# LCs for SED project 2013
     757table="AnalysisResultsRunCutsLC" # CutsLC
     758# Mrk 421
     759source=1
     760# nightly
     761bin=-1
     762nightmin=20120901
     763nightmax=20130516
     764name="Mrk421_nightly_SED-project-2013"
     765get_results
     766
     767
     768
     769# LCs for XMM proposal
     770table="AnalysisResultsRunCutsLC" # CutsLC
     771# Mrk 421
     772source=1
     773# nightly
     774bin=-1
     775nightmin=20111115
     776nightmax=20201231
     777name="Mrk421_nightly_XMM-proposal-2019"
     778get_results
    660779# Mrk 501
    661780source=2
    662 name="Mrk501_all_P"
    663 get_results
    664 
    665 #exit
     781name="Mrk501_nightly_XMM-proposal-2019"
     782get_results
     783# 1959
     784source=7
     785name="1959_nightly_XMM-proposal-2019"
     786get_results
     787
     788
     789
     790# LCs for 421 campaign with Astrosat Jan 2019
     791table="AnalysisResultsRunCutsLC" # CutsLC
     792# Mrk 421
     793source=1
     794# nightly
     795bin=-1
     796nightmin=20190109
     797nightmax=20190113
     798name="Mrk421_nightly_Astrosat-Jan-2019"
     799get_results
     800# 20min
     801bin=20
     802name="Mrk421_20min_Astrosat-Jan-2019"
     803get_results
     804# 30min
     805bin=30
     806name="Mrk421_30min_Astrosat-Jan-2019"
     807get_results
     808# 40min
     809bin=40
     810name="Mrk421_40min_Astrosat-Jan-2019"
     811get_results
     812# 60min
     813bin=60
     814name="Mrk421_60min_Astrosat-Jan-2019"
     815get_results
     816# 80min
     817bin=80
     818name="Mrk421_80min_Astrosat-Jan-2019"
     819get_results
     820# 90min
     821bin=90
     822name="Mrk421_90min_Astrosat-Jan-2019"
     823get_results
     824
     825
     826
     827# LCs for flaring pattern project
     828table="AnalysisResultsRunCutsLC" # CutsLC
     829# Mrk 421
     830source=1
     831# nightly
     832bin=-1
     833nightmin=20111115
     834nightmax=20201231
     835name="Mrk421_nightly_Flaring-pattern-project"
     836get_results
     837# 2d
     838bin=-2
     839nightmin=20111115
     840nightmax=20201231
     841name="Mrk421_2day_Flaring-pattern-project"
     842get_results
     843
     844
     845
     846# LCs for ToO-Trigger project
     847table="AnalysisResultsRunCutsLC" # CutsLC
     848# Mrk 421
     849source=1
     850# nightly for 2019
     851bin=-1
     852nightmin=20190101
     853nightmax=20190630
     854name="Mrk421_nightly_20190101-20190630_ToO-project-2019"
     855get_results
     856# 20 min around flare
     857bin=20
     858nightmin=20190609
     859nightmax=20190612
     860name="Mrk421_20min_20190609-20190612_ToO-project-2019"
     861get_results
     862
     863exit
     864
     865table="AnalysisResultsRunLP" # QLA
     866# Mrk 421
     867source=1
     868# nightly for 2019
     869bin=-1
     870nightmin=20190101
     871nightmax=20190630
     872name="QLA_Mrk421_nightly_20190101-20190630_ToO-project-2019"
     873get_results
     874# 20 min around flare
     875bin=20
     876nightmin=20190609
     877nightmax=20190612
     878name="QLA_Mrk421_20min_20190609-20190612_ToO-project-2019"
     879get_results
     880
     881
     882table="AnalysisResultsRunISDC" # ISDC
     883# Mrk 421
     884source=1
     885# nightly for 2019
     886bin=-1
     887nightmin=20190101
     888nightmax=20190630
     889name="ISDC_Mrk421_nightly_20190101-20190630_ToO-project-2019"
     890get_results
     891# 20 min around flare
     892bin=20
     893nightmin=20190609
     894nightmax=20190612
     895name="ISDC_Mrk421_20min_20190609-20190612_ToO-project-2019"
     896get_results
     897
     898
     899exit
    666900
    667901# LC for ICRC
     
    692926
    693927
    694 #exit
     928exit
    695929
    696930# Mrk 501
Note: See TracChangeset for help on using the changeset viewer.