discovery.sh

Go to the documentation of this file.

#!/bin/env zsh
#
# \author mdejong
#
version=1.0
script=${0##*/}
 
source $JPP_DIR/setenv.sh $JPP_DIR >& /dev/null
 
set_variable   DEBUG       2
set_variable:  NUISANCE    ASTRONOMY_NUISANCE     fixed fixed
set_variable:  SNR         ASTRONOMY_SNR          0.0
set_variable:  M_SIZE      ASTRONOMY_M_SIZE       10000
set_variable:  FORMAT      GRAPHICS_FORMAT        gif
set_variable+  BATCH       GRAPHICS_BATCH         -B
 
if do_usage $*; then
    usage "$script <signal histogram> <background histogram> <number of tests> <probability>"\
          "\nThe histograms correspond to \"<file name>:<histogram name>\"."
fi
 
if (( $# != 4 )); then
    fatal "Wrong number of arguments."
fi
 
set_variable  HS          $argv[1]
set_variable  HB          $argv[2]
let  "NUMBER_OF_TESTS  =  $argv[3]"
let  "PROBABILITY      =  $argv[4]"
 
if (( $PROBABILITY <= 0.0 || $PROBABILITY >= 1.0 )); then
    fatal "Invalid probability $PROBABILITY"
fi
 
if (( $NUMBER_OF_TESTS * $PROBABILITY <= 1.0 )); then
    warning "Insufficient number of tests $NUMBER_OF_TESTS given probability $PROBABILITY."
fi
 
timer_start
 
# H0 pseudo experiments
 
let "N = 0"                  # zero signal
 
H0=$N
 
echo -n "Running H0 pseudo experiments... "
 
$JPP_DIR/examples/JAstronomy/JPseudoExperiment \
    -o benchmark\[$H0\].root \
    -E "$HS $HB"             \
    -N "$NUISANCE"           \
    -n $NUMBER_OF_TESTS      \
    -s $N                    \
    -x "10000 -10 100"       \
    -P $PROBABILITY          \
    -R $SNR                  \
    -M $M_SIZE               \
    -S 0                     \
    -d $DEBUG                >& benchmark.$H0.log
 
echo OK
cat benchmark.$H0.log
 
awk '/Minimal likelihood ratio:/ { print $(NF - 1) }' benchmark.$H0.log | read LIKELIHOOD_RATIO
 
printf "Minimal likelihood ratio: %9.5f\n" $LIKELIHOOD_RATIO
 
 
# H1 pseudo experiments
 
let "NUMBER_OF_TESTS  =  100000"
let "QUANTILE         =  0.5"
 
$JPP_DIR/examples/JAstronomy/JPseudoExperiment \
    -o /dev/null             \
    -E "$HS $HB"             \
    -N "$NUISANCE"           \
    -n $NUMBER_OF_TESTS      \
    -s 1.0                   \
    -Q $LIKELIHOOD_RATIO     \
    -P $QUANTILE             \
    -R $SNR                  \
    -M $M_SIZE               \
    -S 0                     \
    -d $DEBUG                >& benchmark.log
 
awk '/Signal strength:/ { print $NF }' benchmark.log | read X
 
printf "signal: %9.5f\n" $X
 
printf "%1.1f" $X | read H1
 
$JPP_DIR/examples/JAstronomy/JPseudoExperiment \
    -o benchmark\[$H1\].root \
    -E "$HS $HB"             \
    -N "$NUISANCE"           \
    -n $NUMBER_OF_TESTS      \
    -s $X                    \
    -x "1000 -10 100"        \
    -R $SNR                  \
    -M $M_SIZE               \
    -S 0                     \
    -d $DEBUG                >& benchmark.log
 
timer_stop
timer_print
 
 
# check
 
if (( 1 )); then
 
    typeset -A HYPOTHESES
 
    HYPOTHESES[H0]=$PROBABILITY 
    HYPOTHESES[H1]=$QUANTILE
 
    for KEY VALUE in ${(kv)HYPOTHESES}; do
 
        JZoom1D \
            -f benchmark\[${(P)KEY}\].root:hl    \
            -x "$LIKELIHOOD_RATIO 1.0e99"        \
            -o zoom.root
 
        JPrintResult -f benchmark\[${(P)KEY}\].root:hl -F GetSumOfWeights  |  read Y0
        JPrintResult -f zoom.root:hl                   -F GetSumOfWeights  |  read Y1
 
        printf "Hypothesis %s - check probability: %10.3e == %10.3e\n"  $KEY  $VALUE  $(($Y1 / $Y0))
 
        rm -f zoom.root
    done
fi
 
 
# plot
 
let "X1 = $LIKELIHOOD_RATIO"
let "X2 = $LIKELIHOOD_RATIO"
let "Y1 = 1.0"
let "Y2 = $(JPrintResult -f benchmark\[$N\].root:hl -F GetMaximum)"
 
JLine -p "$X1 $Y1 $X2 $Y2" -@ "style = 2" -o line.root
 
JPlot1D \
    -f benchmark\[$H0\].root:hl    \
    -f benchmark\[$H1\].root:hl    \
    -f line.root:\.\*              \
    -Y                             \
    -> "likelihood ratio"         \
    -\^ "number of experiments"    \
    -T ""                          \
    -O HIST                        \
    -L TR                          \
    -o discovery.$FORMAT           $BATCH