plot-waveforms.sh

Go to the documentation of this file.

#!/bin/zsh
script=${0##*/}

if [ -z $JPP_DIR ]; then
    echo "Variable JPP_DIR undefined."
    exit
fi

source $JPP_DIR/setenv.sh $JPP_DIR  >&  /dev/null

set_variable   DEBUG                0
set_variable   WORKDIR              ${TMPDIR:-/tmp}/.$$
set_variable   FREQUENCY_HZ         $((25.0e6 / 128))
set_variable:  OPTION      GRAPHICS_OPTION
set_variable:  FORMAT      GRAPHICS_FORMAT        gif
set_variable+  BATCH       GRAPHICS_BATCH         -B

if do_usage $*; then
    usage "$script (<detector identifier>)+"\
          "\nAuxiliary script to plot frequency spectra of waveforms for given detector identifier(s)."
fi

if (( $# == 0 )); then
    fatal "Wrong number of arguments."
fi

source JAcousticsToolkit.sh         >&  /dev/null

set_array      DETECTORS            $argv[1,-1]
set_variable   ACOUSTICS_WORKDIR    $WORKDIR
set_variable   FORMULA              "sin([0] + 6.28 * ([1] + [2]*x) * x)"
set_variable   DY                   1.0e-8


mkdir   $WORKDIR

for DETECTOR in $DETECTORS[*]; do

    (rm -f $WORKDIR/*.txt) >& /dev/null

    if [[ $DETECTOR == [0-9][0-9]* ]]; then
        set_variable  ID        $DETECTOR
        set_variable  DETECTOR  `getDetector -D $ID`
    else
        set_variable  ID        `getDetector -D $DETECTOR`
    fi

    JAcoustics.sh $ID >& /dev/null

    typeset -A WAVEFORMS

    WAVEFORMS=()

    get_waveforms  $WORKDIR/waveform.txt  WAVEFORMS

    for KEY VALUE in ${(kv)WAVEFORMS}; do

        if (( $KEY > 0 && $VALUE > 0 )); then

            wget \
                --load-cookies=$HOME/.km3netdb_cookie                              \
                --secure-protocol=TLSv1_2                                          \
                "https://km3netdbweb.in2p3.fr/xmlds/acwavfrm/s?waveformid=${KEY}"  \
                -O $WORKDIR/wv.xml                                                 >& /dev/null

            xmllint $WORKDIR/wv.xml                                                                 |\
            awk '/<Description>/,/<\/Description>/       { print gensub("<[^>]*>", "", "g", $0) }'  |  read DESCRIPTION

            xmllint $WORKDIR/wv.xml                                                                 |\
            awk '/<EncodedSamples>/,/<\/EncodedSamples>/ { print gensub("<[^>]*>", "", "g", $0) }'  |\
            tail -n +2                                                                              >  $WORKDIR/wv.txt


            JFFT \
                -f $WORKDIR/wv.txt                                 \
                -o $WORKDIR/fft.txt                                \
                -B $((1.0e+3 / $FREQUENCY_HZ))

            JGraph \
                -f $WORKDIR/fft.txt                                \
                -o $WORKDIR/${DETECTOR}\[${VALUE}\].fft.root

            JPrintQuantiles \
                -f $WORKDIR/${DETECTOR}\[${VALUE}\].fft.root:\.\*  \
                -Q "0.3 0.7"                                       |  read -A RESULT
                
            printf "%s %2d %7.3f %7.3f [kHz] \"%s\"\n" $DETECTOR $VALUE $((0.5 * ($RESULT[1] + $RESULT[2]))) $(($RESULT[2] - $RESULT[1])) $DESCRIPTION

            Y=(`cat $WORKDIR/wv.txt`)

            rm -f $WORKDIR/graph.txt

            for (( i = 1; $i <= ${#Y}; ++i )); do

                X=$((($i - 1) / $FREQUENCY_HZ))

                echo  $X  $Y[$i]  $DY  >>  $WORKDIR/graph.txt
            done

            JGraph \
                -f $WORKDIR/graph.txt                                       \
                -o $WORKDIR/graph.root

            # refine start values

            JFit \
                -f $WORKDIR/graph.root:\.\*                                 \
                -o $WORKDIR/fit\[${KEY}\].root                              \
                -x "0 0.1e-3"                                               \
                -F "$FORMULA"                                               \
                -= "p0 = 0.0;"                                              \
                -@ "p1 = $((1.0e3 *  $RESULT[1]));"                         \
                -= "p2 = 0.0;"                                              \
                -d $DEBUG

            JPrintFit -f $WORKDIR/fit\[${KEY}\].root:\.\* -@ p0  |  read -A P0
            JPrintFit -f $WORKDIR/fit\[${KEY}\].root:\.\* -@ p1  |  read -A P1
            JPrintFit -f $WORKDIR/fit\[${KEY}\].root:\.\* -@ p2  |  read -A P2

            JFit \
                -f $WORKDIR/graph.root:\.\*                                \
                -o $WORKDIR/fit\[${KEY}\].root                             \
                -x "0 0.5e-3"                                              \
                -F "$FORMULA"                                              \
                -= "p0 = $P0[1];"                                          \
                -@ "p1 = $P1[1];"                                          \
                -@ "p2 = $P2[1];"                                          \
                -d $DEBUG

            JPrintFit -f $WORKDIR/fit\[${KEY}\].root:\.\* -@ p0  |  read -A P0
            JPrintFit -f $WORKDIR/fit\[${KEY}\].root:\.\* -@ p1  |  read -A P1
            JPrintFit -f $WORKDIR/fit\[${KEY}\].root:\.\* -@ p2  |  read -A P2
            
            JFit \
                -f $WORKDIR/graph.root:\.\*                                \
                -o $WORKDIR/fit\[${KEY}\].root                             \
                -F "$FORMULA"                                              \
                -= "p0 = $P0[1];"                                          \
                -@ "p1 = $P1[1];"                                          \
                -@ "p2 = $P2[1];"                                          \
                -d $DEBUG

            JPrintFit -f $WORKDIR/fit\[${KEY}\].root:\.\* -@ p0                     |  read -A P0
            JPrintFit -f $WORKDIR/fit\[${KEY}\].root:\.\* -@ p1                     |  read -A P1
            JPrintFit -f $WORKDIR/fit\[${KEY}\].root:\.\* -@ p2                     |  read -A P2
            JPrintFit -f $WORKDIR/fit\[${KEY}\].root:\.\* -F "GetChisquare/GetNDF"  |  read CHI2
            
            printf "%s %2d %7.3f [kHz] %9.3f [kHz/s] %12.3e\n" $DETECTOR $KEY $((1.0e-3 * $P1[1])) $((1.0e-3 * $P2[1])) $CHI2

            JPlot1D \
                -w 1200x600                                                 \
                -f $WORKDIR/fit\[${KEY}\].root:\.\*                         \
                -> "time [s]"                                              \
                -\^ "amplitude [a.u.]"                                      \
                -N "X 505"                                                  \
                -L TR                                                       \
                -T "$DETECTOR ($ID)"                                        \
                -o ${DETECTOR}.${KEY}.$FORMAT                               $BATCH
        fi
    done

    INPUT_FILES=(`ls $WORKDIR/${DETECTOR}\[*\].fft.root`)

    JPlot1D \
        -f${^INPUT_FILES}:\.\*                                     \
        -x "10 50"                                                 \
        -> "frequency [kHz]"                                      \
        -\^ "power [a.u.]"                                         \
        $OPTION                                                    \
        -L TR                                                      \
        -T "$DETECTOR ($ID)"                                       \
        -o ${DETECTOR}.$FORMAT                                     $BATCH
done
    
rm -rf  $WORKDIR