hydrophone-t0:fit.sh

Go to the documentation of this file.

#!/usr/bin/env zsh
script=${0##*/}
 
if [ -z $JPP_DIR ]; then
    echo "Variable JPP_DIR undefined."
    exit
fi
 
source $JPP_DIR/setenv.sh $JPP_DIR >& /dev/null
 
zmodload zsh/mathfunc
 
set_variable:  DEBUG       ACOUSTICS_DEBUG        2
set_variable:  WORKDIR     ACOUSTICS_WORKDIR      ./
set_variable:  XMIN        HYDROPHONE_XMIN        -1.0e-3
set_variable:  XMAX        HYDROPHONE_XMAX        +1.0e-3
set_variable:  FORMAT      GRAPHICS_FORMAT        gif
set_variable+  BATCH       GRAPHICS_BATCH         -B
 
if do_usage $*; then
    usage "$script <detector file> <string identifier> [<detector file>]"
fi
 
if (( $# == 0 || $# > 3 )); then
    fatal "Wrong number of arguments."
fi
 
set_variable  DETECTOR            $argv[1]
set_variable  STRING              $argv[2]
set_array     QUANTILES           0.333 0.666
set_variable  FORMULA             "[0] * exp(-0.5 * (x-[1])*(x-[1]) / ([2]*[2]))"
set_variable  MODULE              `getModule -a $DETECTOR -L "$STRING 0"`
 
source JAcousticsToolkit.sh
 
typeset -A TRIPODS
 
get_tripods $WORKDIR/tripod.txt TRIPODS
 
XMEAN=(0.0 0.0)
STDEV=(0.0 0.0)
TOTAL=0.0
 
for TRIPOD in ${(k)TRIPODS[*]}; do
    
    JPrintQuantiles \
        -f "hydrophone.root:${MODULE}.${TRIPOD}"     \
        -Q "$QUANTILES[*]"                           \
        -d 0                                         | read -A RESULT
 
    if (( ${#RESULT} == ${#QUANTILES} )); then
        
        P1=$((0.5*($RESULT[-1] + $RESULT[1])))           # average
        P2=$((0.5*($RESULT[-1] - $RESULT[1])))           # width
 
        JFit \
            -f "hydrophone.root:${MODULE}.${TRIPOD}"     \
            -o ${TMPDIR:-/tmp}/fit.root                  \
            -F "$FORMULA"                                \
            -@"p0 = GetMaximum * 0.8"                    \
            -@"p1 = $P1"                                 \
            -@"p2 = $P2"                                 \
            -d 0                                         >& /dev/null
 
        JPrintFit -f  ${TMPDIR:-/tmp}/fit.root:\.\* -@p1  |  read -A P1
        JPrintFit -f  ${TMPDIR:-/tmp}/fit.root:\.\* -@p2  |  read -A P2
 
        VALUE=$((1.0e+6 * $P1[1]))                              # [us]
        TOTAL=$(($TOTAL + 1.0))
 
        if (( $TOTAL == 1.0 )); then
            XMEAN[1]=$VALUE
            STDEV[1]=0.0
        else
            XMEAN[2]=$(($XMEAN[1] + ($VALUE - $XMEAN[1]) / $TOTAL))
            STDEV[2]=$(($STDEV[1] + ($VALUE - $XMEAN[1]) * ($VALUE - $XMEAN[2])))
            XMEAN[1]=$XMEAN[2]
            STDEV[1]=$STDEV[2]
        fi
 
        printf "string %4d tripod %2d %5.0f %3.0f us\n" $STRING $TRIPOD $(($P1[1] * 1.0e6)) $(($P2[1] * 1.0e6))
    fi
done
 
if (( $TOTAL > 1.0 )); then
    printf "string %4d mean  RMS %5.0f %3.0f us\n" $STRING $XMEAN[1] $((sqrt($TOTAL * $STDEV[1] / ($TOTAL * ($TOTAL - 1.0)))))
fi
 
typeset -Z 4 STRING
 
JOpera1D \
    -f hydrophone.root:$MODULE                   \
    -o ${TMPDIR:-/tmp}/H\[${STRING}\].root       \
    -u Add -d 0
 
JPrintQuantiles \
    -f ${TMPDIR:-/tmp}/H\[${STRING}\].root:Add   \
    -Q "$QUANTILES[*]"                           \
    -d 0                                         | read -A RESULT
 
P1=$((0.5*($RESULT[-1] + $RESULT[1])))           # average
P2=$((0.5*($RESULT[-1] - $RESULT[1])))           # width
 
JFit \
    -f ${TMPDIR:-/tmp}/H\[${STRING}\].root:Add   \
    -o ${TMPDIR:-/tmp}/F\[${STRING}\].root       \
    -F "$FORMULA"                                \
    -@"p0 = GetMaximum * 0.8"                    \
    -@"p1 = $P1"                                 \
    -@"p2 = $P2"                                 \
    -d $DEBUG
 
JPrintFit -f  ${TMPDIR:-/tmp}/F\[${STRING}\].root:Add -@p1  |  read -A P1
JPrintFit -f  ${TMPDIR:-/tmp}/F\[${STRING}\].root:Add -@p2  |  read -A P2
 
printf "string %4d t0    %6.0f +/- %6.0f us\n" $STRING $(($P1[1] * 1.0e6)) $(($P1[2] * 1.0e6))
printf "string %4d sigma %6.0f +/- %6.0f us\n" $STRING $(($P2[1] * 1.0e6)) $(($P2[2] * 1.0e6))
 
if (( ${#} == 3 )); then
 
    set_variable  DETECTOR            $argv[3]
 
    set_array P1  `JPrintFit -f  ${TMPDIR:-/tmp}/F\[${STRING}\].root:Add -@p1`
 
    let "TIME_OFFSET_NS = $P1[1] * 1.0e9"
 
    printf "Add time offset %12.6d to hydrophone.\n" $TIME_OFFSET_NS
    
    JEditDetector \
        -a $DETECTOR                         \
        -o $DETECTOR                         \
        -M "$MODULE SUB $TIME_OFFSET_NS"     \
        -d $DEBUG
fi
 
JPlot1D \
    -f ${TMPDIR:-/tmp}/F\[${STRING}\].root:Add   \
    -x "$XMIN $XMAX"                             \
    -> "#Deltat [s]"                            \
    -\^ "number of events [a.u.]"                \
    -N "X 505"                                   \
    -T "$STRING"                                 \
    -o hydrophonet0_${STRING}.$FORMAT            $BATCH