examples/scripts/gradient.sh

Go to the documentation of this file.

#!/usr/bin/env zsh
script=${0##*/}
 
source $JPP_DIR/setenv.sh $JPP_DIR >& /dev/null
source $JPP_BIN/gradient.sh
 
set_variable:  DEBUG       SCRIPTS_DEBUG           3
 
if do_usage $*; then
    usage "$script [number of tests]"
fi
 
if (( $# == 1 )); then
    set_variable NUMBER_OF_TESTS      $1
else
    set_variable NUMBER_OF_TESTS      0
fi
 
#
# User function.
#
# \param 1     x
# \param 2     mean
# \param 3     sigma
# \param 4     signal
# \param 5     background
#
function gauss()
{
    let "PI = acos(-1.0)"
    let "u  = ($1 - $2) / $3"
 
    echo $(( $4 * exp(-0.5*u*u) / (sqrt(2.0*$PI) * $3)  +  $5 ))
}
 
 
#
# User data.
#
typeset -A    DATA                # (x,y) values
 
 
#
# Generate random data according true parameters.
#
function generate
{
    DATA=()
 
    for (( X = -5.0; $X <= +5.0; X += 0.5 )); do
 
        let "Y  = `gauss $X $TRUE_MEAN $TRUE_SIGMA $TRUE_SIGNAL $TRUE_BACKGROUND`"
        
        DATA[$X]=`$JPP_DIR/examples/JMath/getPoisson -e $Y`
    done
}
 
 
#
# Get chi2.
#
# The function is evaluated using parameters MEAN, SIGMA, SIGNAL and BACKGROUND.
#
# \param  1       variable containing chi2 value on return
# \param  2       option
#
function getChi2()
{
    let "SUM = 0.0"
 
    for X Y in ${(kv)DATA}; do
 
        let "F  = `gauss $X $MEAN $SIGMA $SIGNAL $BACKGROUND`"
        let "DY = ($F - $Y) * ($F - $Y)"
 
        if (( $Y > 0 )); then
            let "DY /= $Y"
        else
            let "DY  = $F"
        fi
 
        let "SUM += $DY"
 
    done
 
    eval $1=$SUM
}
 
 
# True values.
 
let "TRUE_MEAN        =    0.0"
let "TRUE_SIGMA       =    1.0"
let "TRUE_SIGNAL      = 1000.0"
let "TRUE_BACKGROUND  =   10.0"
 
 
#
# Print parameters.
#
function typeout()
{
    for KEY in MEAN SIGMA SIGNAL BACKGROUND; do
 
        set_variable P1 $KEY
        set_variable P2 TRUE_$KEY
        
        printf "%-12s  %12.3f (fit) %12.3f (true)\n"  $KEY  ${(P)P1}  ${(P)P2}
    done
}
 
 
if (( $NUMBER_OF_TESTS == 0 )); then
 
    DATA=(
        -3.00000     16.00000
        -2.50000     22.00000
        -2.00000     70.00000
        -1.50000    152.00000
        -1.00000    261.00000
        -0.50000    337.00000
         0.00000    378.00000
         0.50000    360.00000
         1.00000    253.00000
         1.50000    129.00000
         2.00000     72.00000
         2.50000     22.00000
         3.00000     11.00000)
 
    # Fit parameters.
 
    PARAMETERS=()
    
    PARAMETERS['let "MEAN        +=  %"']=1.0e-2
    PARAMETERS['let "SIGMA       +=  %"']=1.0e-2
    PARAMETERS['let "SIGNAL      +=  %"']=5.0e-0
    PARAMETERS['let "BACKGROUND  +=  %"']=5.0e-1
 
    # Start values.
 
    let "MEAN        =    0.5"
    let "SIGMA       =    0.5"
    let "SIGNAL      =  700.0"
    let "BACKGROUND  =    0.0"
 
    gradient 0  Y1
 
    printf "chi2/NDF      %12.5f (fit) (gradient)\n" $(($Y1 / (${#DATA} - ${#PARAMETERS})))
 
    simplex  Y1
 
    printf "chi2/NDF      %12.5f (fit) (simplex) \n" $(($Y1 / (${#DATA} - ${#PARAMETERS})))
 
    typeout
 
    let "MEAN        = $TRUE_MEAN"
    let "SIGMA       = $TRUE_SIGMA"
    let "SIGNAL      = $TRUE_SIGNAL"
    let "BACKGROUND  = $TRUE_BACKGROUND"
 
    getChi2  Y2
 
    printf "chi2/NDF      %12.5f (fit) %12.5f (true)\n" $(($Y1 / (${#DATA} - ${#PARAMETERS})))  $(($Y2 / ${#DATA}))
 
else
 
    set_variable:  DEBUG       SCRIPTS_DEBUG           1
 
    let "CHI2_FIT  = 0.0"
    let "CHI2_TRUE = 0.0"
 
    for (( m = 0; $m != $NUMBER_OF_TESTS; ++m )); do
 
        generate
 
        # Fit parameters.
 
        PARAMETERS=()
    
        PARAMETERS['let "MEAN        +=  %"']=1.0e-2
        PARAMETERS['let "SIGMA       +=  %"']=1.0e-2
        PARAMETERS['let "SIGNAL      +=  %"']=5.0e-0
        PARAMETERS['let "BACKGROUND  +=  %"']=5.0e-1
 
        # Start values.
 
        let "MEAN        =    0.5"
        let "SIGMA       =    0.5"
        let "SIGNAL      =  700.0"
        let "BACKGROUND  =    0.0"
 
        gradient 0  Y1
 
        let "MEAN        = $TRUE_MEAN"
        let "SIGMA       = $TRUE_SIGMA"
        let "SIGNAL      = $TRUE_SIGNAL"
        let "BACKGROUND  = $TRUE_BACKGROUND"
        
        getChi2  Y2
        
        printf "chi2/NDF  %3d %12.5f (fit) %12.5f (true)\n" $m  $(($Y1 / (${#DATA} - ${#PARAMETERS})))  $(($Y2 / ${#DATA}))
 
        let "CHI2_FIT  += $Y1 / (${#DATA} - ${#PARAMETERS})"
        let "CHI2_TRUE += $Y2 /  ${#DATA}"
    done
 
    printf "<>            %12.5f (fit) %12.5f (true)\n"  $(($CHI2_FIT / $NUMBER_OF_TESTS))  $(($CHI2_TRUE / $NUMBER_OF_TESTS))
fi