observation_mset.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
 
zmodload zsh/mathfunc
 
set_variable   DEBUG       2
set_variable:  SNR         ASTRONOMY_SNR          0.2
set_variable:  M_SIZE      ASTRONOMY_M_SIZE       10000
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 <config file name for data> <config file name for pseudosets>"\
          "\nThe histograms correspond to \"<file name>:<histogram name>\"."
fi
 
if (( $# != 2 )); then
    fatal "Wrong number of arguments."
fi
 
set_variable  CONFIG_DATA       $argv[1]
set_variable  CONFIG_PSEUDO     $argv[2]
 
timer_start
 
let "NUMBER_OF_TESTS  =  100000"
 
$JPP_DIR/examples/JAstronomy/JRealExperiment \
        -E $CONFIG_DATA             \
        -R $SNR                     \
        -n $NUMBER_OF_TESTS         \
        -d $DEBUG                   >& benchmark.data.log
 
awk '/result:/      { print $(NF - 1) }' benchmark.data.log | read LIKELIHOOD_RATIO
awk '/result:/      { print $(NF) }'     benchmark.data.log | read SIGNAL
awk '/upper limit:/ { print $(NF) }'     benchmark.data.log | read UPPER_LIMIT
 
SIGMA_UPPER_LIMIT=$(( $UPPER_LIMIT * sqrt(1.0 / ($NUMBER_OF_TESTS*1.0)*0.9/0.1) ))
 
#test p-value calculation: for a best fitted signal build pseudoexperiments with such a signal and find fraction of pex with TS > TS_observed
$JPP_DIR/examples/JAstronomy/JGen2 \
    -o /dev/null             \
    -E $CONFIG_PSEUDO        \
    -n 10000                 \
    -s $SIGNAL               \
    -Q $LIKELIHOOD_RATIO     \
    -R $SNR                  \
    -M $M_SIZE               \
    -d $DEBUG                >& benchmark.$SIGNAL.log
 
awk '/Maximal probability:/ { print $(NF - 1) }' benchmark.$SIGNAL.log | read PCHECK
 
if (( $(echo "$LIKELIHOOD_RATIO <= 0" |bc -l) )); then
        printf "Observed likelihood ratio: %1.2e signal\t%1.2e pcheck\t%1.2e\n\t\t\t\t\t\t" $LIKELIHOOD_RATIO $SIGNAL $PCHECK
fi
 
H0=0
 
let "NUMBER_OF_TESTS  =  100000"
 
if (( $LIKELIHOOD_RATIO > 0 )); then
 
        for (( ; 1 ; )); do              # p-vale
                $JPP_DIR/examples/JAstronomy/JGen2 \
                    -o /dev/null             \
                    -E $CONFIG_PSEUDO        \
                    -n $NUMBER_OF_TESTS      \
                    -s $H0                   \
                    -Q $LIKELIHOOD_RATIO     \
                    -R $SNR                  \
                    -M $M_SIZE               \
                    -d $DEBUG                >& benchmark.$H0.log
 
            awk '/Maximal probability:/ { print $(NF - 1) }' benchmark.$H0.log | read P
 
        if (( $P*($NUMBER_OF_TESTS*1.0) > 1.0e3 )); then 
                    break
            else
                    let "NUMBER_OF_TESTS *= 10"
            fi
        done
 
        STDERR=$(( sqrt($P/($NUMBER_OF_TESTS*1.0)) ))
 
        printf "Observed likelihood ratio: %1.2e signal\t%1.2e pcheck\t%1.2e\t -> P = %1.2e+-%1.2e, \t n=%1.2e\tUL = %1.2e+-%1.2e \n" $LIKELIHOOD_RATIO $SIGNAL $PCHECK $P $STDERR $NUMBER_OF_TESTS $UPPER_LIMIT $SIGMA_UPPER_LIMIT
fi