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
 
set_variable   DEBUG       2
set_variable:  NUISANCE    ASTRONOMY_NUISANCE     fixed fixed
set_variable:  FORMAT      GRAPHICS_FORMAT        gif
set_variable+  BATCH       GRAPHICS_BATCH         -B
set_variable:  SNR         ASTRONOMY_SNR          0.2
set_variable:  M_SIZE      ASTRONOMY_M_SIZE       10000
 
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
 
 
$JPP_DIR/examples/JAstronomy/JRealExperiment \
        -E "`cat $CONFIG_DATA`"     \
        -R $SNR                     \
        -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
 
#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 benchmark\[$H0\].root \
    -E $CONFIG_PSEUDO        \
    -n 10000                 \
    -s $SIGNAL               \
    -x "10000 -10 100"       \
    -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 (( $(echo "$LIKELIHOOD_RATIO > 0" |bc -l) )); then
        for (( ; 1 ; )); do              # p-vale
                $JPP_DIR/examples/JAstronomy/JGen2 \
                    -o benchmark\[$H0\].root \
                    -E $CONFIG_PSEUDO        \
                    -n $NUMBER_OF_TESTS      \
                    -s 0                     \
                    -x "10000 -10 100"       \
                    -Q $LIKELIHOOD_RATIO     \
                    -R $SNR                  \
                    -M $M_SIZE               \
                    -d $DEBUG                >& benchmark.$H0.log
 
            awk '/Maximal probability:/ { print $(NF - 1) }' benchmark.$H0.log | read P
 
            NTAIL_OK=$(awk -v NUMBER_OF_TESTS="$NUMBER_OF_TESTS" -v P="$P" 'BEGIN {print (P*NUMBER_OF_TESTS > 1e3 ? 1 : 0)}')
        
            if (( $NTAIL_OK == 1 )); then
                break
            else
                let "NUMBER_OF_TESTS *= 10"
            fi
        done
        STDERR=$(awk -v NUMBER_OF_TESTS="$NUMBER_OF_TESTS" -v P="$P" 'BEGIN {print (sqrt(P/NUMBER_OF_TESTS))}')
 
        printf "Observed likelihood ratio: %1.2e signal\t%1.2e pcheck\t%1.2e\t -> P = %1.2e+-%1.2e, \t n=%1.2e\n" $LIKELIHOOD_RATIO $SIGNAL $PCHECK $P $STDERR $NUMBER_OF_TESTS
fi