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:  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
 
# pseudo experiments
 
let "N = 0"                  # zero signal
 
$JPP_DIR/examples/JAstronomy/JPseudoExperiment \
    -o benchmark\[$N\].root  \
    -E "$HS $HB"             \
    -n $NUMBER_OF_TESTS      \
    -s $N                    \
    -x "10000 -10 100"       \
    -N "$NUISANCE"           \
    -P $PROBABILITY          \
    -d $DEBUG --!            >& benchmark.$N.log
 
cat benchmark.$N.log
 
awk '/Minimal likelihood ratio:/ { print $(NF - 1) }' benchmark.$N.log | read LIKELIHOOD_RATIO
 
printf "Minimal likelihood ratio: %9.5f\n" $LIKELIHOOD_RATIO
 
 
# check
 
JZoom1D \
    -f benchmark\[$N\].root:hl      \
    -x "$LIKELIHOOD_RATIO 1.0e99"   \
    -o zoom\[$N\].root
 
JPrintResult -f benchmark\[$N\].root:hl -F GetSumOfWeights  |  read Y0
JPrintResult -f zoom\[$N\].root:hl      -F GetSumOfWeights  |  read Y1
 
printf "Check of probability: %12.3e  %12.3e\n" $PROBABILITY  $(($Y1 / $Y0))
 
 
# binary search
 
let "NUMBER_OF_TESTS  =  100000"
 
let "XMIN      =   1.0"
let "XMAX      =   2.0"
let "EPS       =   1.0e-2"
let "QUANTILE  =   0.5"
 
for (( ; 1 ; )); do              # lower limit
 
    $JPP_DIR/examples/JAstronomy/JPseudoExperiment \
        -o /dev/null             \
        -E "$HS $HB"             \
        -n $NUMBER_OF_TESTS      \
        -s $XMIN                 \
        -Q $LIKELIHOOD_RATIO     \
        -d $DEBUG                >& benchmark.log
 
    awk '/Maximal probability:/ { print $(NF - 1) }' benchmark.log | read PROBABILITY
 
    printf "lower limit: %9.5f -> P = %9.5f\n" $XMIN $PROBABILITY
 
    if (( $PROBABILITY > $QUANTILE )); then
        let "XMAX  = $XMIN"
        let "XMIN *= 0.5"
    else
        break
    fi
done
 
for (( ; 1 ; )); do              # upper limit
 
    $JPP_DIR/examples/JAstronomy/JPseudoExperiment \
        -o /dev/null             \
        -E "$HS $HB"             \
        -n $NUMBER_OF_TESTS      \
        -s $XMAX                 \
        -Q $LIKELIHOOD_RATIO     \
        -d $DEBUG                >& benchmark.log
 
    awk '/Maximal probability:/ { print $(NF - 1) }' benchmark.log | read PROBABILITY
 
    printf "upper limit: %9.5f -> P = %9.5f\n" $XMAX $PROBABILITY
 
    if (( $PROBABILITY < $QUANTILE )); then
        let "XMIN  = $XMAX"
        let "XMAX *= 2.0"
    else
        break
    fi
done
 
for (( ; $XMAX - $XMIN > $EPS; )); do
 
    let "X = 0.5 * ($XMIN + $XMAX)"
 
    $JPP_DIR/examples/JAstronomy/JPseudoExperiment \
        -o /dev/null             \
        -E "$HS $HB"             \
        -n $NUMBER_OF_TESTS      \
        -s $X                    \
        -Q $LIKELIHOOD_RATIO     \
        -d $DEBUG                >& benchmark.log
 
    #cat benchmark.log
 
    awk '/Maximal probability:/ { print $(NF - 1) }' benchmark.log | read PROBABILITY
 
    printf "search: %9.5f <= %9.5f <= %9.5f -> P = %9.5f\n" $XMIN $X $XMAX $PROBABILITY
 
    if (( $PROBABILITY < $QUANTILE )); then
        let "XMIN = $X"
    else
        let "XMAX = $X"
    fi
done
 
let "X = 0.5 * ($XMIN + $XMAX)"
 
$JPP_DIR/examples/JAstronomy/JPseudoExperiment \
    -o benchmark.root        \
    -E "$HS $HB"             \
    -n $NUMBER_OF_TESTS      \
    -s $X                    \
    -x "1000 -10 100"        \
    -S 0                     \
    -d $DEBUG                >& benchmark.log
 
printf "signal: %9.5f\n" $X
 
timer_stop
timer_print