acoustics-fit-toolkit.sh

Go to the documentation of this file.

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

# -----------------------------------------------------------------------------------------------------
#
# Toolkit for the determination of the fixed parameters in the system, a.k.a. "pre-calibration".
# The list of parameters includes (x,y,z) positions of strings and tripods, stretching of strings,
# z-positions of the optical modules, orientations of the strings and possibly z-positions of anchors.
#
# The value of a given parameter is optimised by minimising the average chi2/NDF
# of a series of global fits applied to a predefined set of acoustic events.
# For each parameter, there exists a designated function.
#
# There are also steering functions to consistently determine the optimal values for a set of parameters.
# In these, the current accuracy of the parameter values is taken into account.
# These are therefore referred to as "stages".
# The pre-calibration thus consists of a sequence of stages.
# For each detector, a complete pre-calibration procedure is implemented in a designated script.
#
# Mandatory input variables:
#
# DETECTOR              # detector file
# INPUT_FILES           # list of input files (output of JAcousticsEventBuilder[.sh])
#
# Optional input variables:
#
# HYDROPHONE            # if defined, use hydrophones in global fit; else do not use hydrophones
#
# In the following:
#
# STRINGS corresponds to list of strings to be fitted; this can be changed with function fixStrings.
# TRIPODS corresponds to list of tripods to be fitted; this can be changed with function fixTripods.
#
# TRANSMITTERS corresponds to list of strings which have an emitter on the anchor (maybe empty).
# For these strings, the orientations and the z-positions of the anchors will separately be fitted.
# This list will also be changed with function fixStrings.
# -----------------------------------------------------------------------------------------------------

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

set_variable   TMPDIR     `pwd`/.$$                     # temporary directory
set_variable:  DEBUG      ACOUSTICS_DEBUG      3        # debug level
set_variable:  JOBS       ACOUSTICS_JOBS       0        # number of parallel threads

if do_usage $*; then
    usage "Toolkit for pre-calibrations of detectors."\
          "\nThis script should be sourced by the steering script for a specific detector."\
          "\nConsult the script '$0' for more details."
fi

if [[ ! $ZSH_EVAL_CONTEXT =~ :file$ ]]; then
    echo "This script should be sourced by the steering script for a specific detector."
    exit
fi

if [[ `ps -o stat= -p $$` == *"+" ]]; then
    echo "This script cannot interactively be used."
    exit
fi

source JAcousticsToolkit.sh

typeset -A     TRIPODS               # tripods
typeset -a     CHI2                  # chi2 values; index [1] corresponds to best value

let "NUMBER_OF_ITERATIONS = 100"     # maximal number of iterations
let "EPSILON              = 5.0E-4"  # maximal distance to minimum of chi2
let "RADIUS_M             = 1.0"     # maximal horizontal distance between T-bar and emitter/hydrophone
let "PRECISION            = 10"      # number of decimals of chi2 evaluation


# -----------------------------------------------------------------------------------------------------
# Initialise.
#
# Install input files in working directory and create temporary directory TMPDIR.
#
# -----------------------------------------------------------------------------------------------------
function initialise()
{
    eval `JPrintDetector -a $DETECTOR -O IDENTIFIER`
    eval `JPrintDetector -a $DETECTOR -O SUMMARY`
    eval `JPrintDetector -a $DETECTOR -O CAN`

    JAcoustics.sh  $DETECTOR_ID

    get_tripods  tripod.txt  TRIPODS

    TRANSMITTERS=(`get_strings_with_transmitter transmitter.txt`)

    if [[ "${HYDROPHONE+YES}" == "YES" ]]; then
        HYDROPHONES=(`get_strings_with_hydrophone hydrophone.txt`)
    else
        HYDROPHONES=()
    fi

    mkdir -p  $TMPDIR

    if [[ ! -f pmt.datx ]]; then
        JConvertDetectorFormat -a $DETECTOR -o pmt.datx
    fi
    JDetachPMTs -a $DETECTOR -o $DETECTOR
}


# -----------------------------------------------------------------------------------------------------
# Backup files in given directory.
#
# \param  1       directory
# -----------------------------------------------------------------------------------------------------
function backup()
{
    debug "Backup intermediate results to $1"

    mkdir -p $1

    cp -p  {tripod.txt,transmitter.txt,hydrophone.txt}  $1

    JAttachPMTs -a $DETECTOR -f pmt.datx -o $1/detector.datx
}


# -----------------------------------------------------------------------------------------------------
# Remove temporary directory and restore input files in working directory.
#
# -----------------------------------------------------------------------------------------------------
function clean()
{
    JAttachPMTs -a $DETECTOR -f pmt.datx -o $DETECTOR

    rm -rf $TMPDIR >& /dev/null

    rm -f  {acoustics_fit_parameters,acoustics_trigger_parameters}.txt

    JAcoustics.sh  $DETECTOR_ID
}


# -----------------------------------------------------------------------------------------------------
# Fix strings.
#
# \param  1-N     identifiers
# -----------------------------------------------------------------------------------------------------
function fixStrings()
{
    typeset -a BUFFER

    BUFFER=(`echo $*`)

    if (( ${#BUFFER} > 0 )); then
        STRINGS=(${STRINGS:|BUFFER})
        TRANSMITTERS=(${TRANSMITTERS:|BUFFER})
        HYDROPHONES=(${HYDROPHONES:|BUFFER})
    fi
}


# -----------------------------------------------------------------------------------------------------
# Fix tripods.
#
# \param  1-N     identifiers
# -----------------------------------------------------------------------------------------------------
function fixTripods()
{
    typeset -a BUFFER

    BUFFER=(`echo $*`)
    
    if (( ${#BUFFER} > 0 )); then
        for ID in $BUFFER[*]; do
            unset "TRIPODS[${ID}]"
        done
    fi
}


# -----------------------------------------------------------------------------------------------------
# Evaluate current chi2.
#
# \param  1       variable containing chi2 value on return
# \param  2       option
# -----------------------------------------------------------------------------------------------------
function getChi2()
{
    if   (( $JOBS == 0 )); then

        rm -f $TMPDIR/katoomba.root >& /dev/null

        JKatoomba \
            -a $DETECTOR                                  \
            -f "$INPUT_FILES[*]"                          \
            -o $TMPDIR/katoomba.root                      \
            -T tripod.txt                                 \
            -Y transmitter.txt                            \
            -V sound_velocity.txt                         \
            -M mechanics.txt                              \
            -@ $TMPDIR/acoustics_fit_parameters.txt       \
            -! disable.txt                                \
            ${HYDROPHONE+-H hydrophone.txt}               \
            -d 0                                          >& $TMPDIR/job.log

        if (( $? != 0 )); then
            fatal "non-zero exit code, see $TMPDIR/job.log"
        fi

        set_variable  RESULT  $($JPP_DIR/examples/JAcoustics/JPrintChi2 -f $TMPDIR/katoomba.root -p $PRECISION)

    else

        if   (( $2 == 0 )); then

            #  step wise improvement of the chi2

            {{ JFremantle \
                   -a $DETECTOR                               \
                   -f "$INPUT_FILES[*]"                       \
                   -T tripod.txt                              \
                   -Y transmitter.txt                         \
                   -V sound_velocity.txt                      \
                   -M mechanics.txt                           \
                   -@ $TMPDIR/acoustics_fit_parameters.txt    \
                   -! disable.txt                             \
                   ${HYDROPHONE+-H hydrophone.txt}            \
                   -N $JOBS                                   \
                   -d 0 3>&1 1>&4 } | read RESULT } 4>&1      >& $TMPDIR/job.log

        elif (( $2 == 1 )); then

            #  evaluation of the chi2 before the determination of the gradient of the chi2

            rm -f $TMPDIR/fremantle.jpp >& /dev/null

            {{ JFremantle \
                   -a $DETECTOR                               \
                   -f "$INPUT_FILES[*]"                       \
                   -o $TMPDIR/fremantle.jpp                   \
                   -T tripod.txt                              \
                   -Y transmitter.txt                         \
                   -V sound_velocity.txt                      \
                   -M mechanics.txt                           \
                   -@ $TMPDIR/acoustics_fit_parameters.txt    \
                   -! disable.txt                             \
                   ${HYDROPHONE+-H hydrophone.txt}            \
                   -N $JOBS                                   \
                   -d 0 3>&1 1>&4 } | read RESULT } 4>&1      >& $TMPDIR/job.log

        elif (( $2 == 2 )); then
            
            #  evaluation of the derivative of the chi2 to each fit parameter

            {{ JPlatypus \
                   -a $DETECTOR                               \
                   -f $TMPDIR/fremantle.jpp                   \
                   -T tripod.txt                              \
                   -Y transmitter.txt                         \
                   -V sound_velocity.txt                      \
                   -M mechanics.txt                           \
                   -@ $TMPDIR/acoustics_fit_parameters.txt    \
                   ${HYDROPHONE+-H hydrophone.txt}            \
                   -N $JOBS                                   \
                   -d 0 3>&1 1>&4 } | read RESULT } 4>&1      >& $TMPDIR/job.log

        else

            fatal "invalid option <$2>"
        fi

        if (( $? != 0 )); then
            fatal "non-zero exit code, see $TMPDIR/job.log"
        fi
    fi

    let "UV = $RESULT"

    if (( $? != 0 )); then
        fatal "non-zero exit code let, value <${RESULT}>."
    fi

    eval $1=$UV
}


# -----------------------------------------------------------------------------------------------------
# Stage 0.
# This stage can be used to determine the (x,y,z) positions of tripods using strings that are fixed.
#
# \param  1       fit option
# \param  2       M-estimator
# \param  3       resolution [s]
# \param  4       standard deviations
# \param  5       step size  [m]
# \param  6       maximum number of extra steps
# -----------------------------------------------------------------------------------------------------
function stage_0()
{
    cat>$TMPDIR/acoustics_fit_parameters.txt<<EOF
`egrep Tmax_s acoustics_fit_parameters.txt`
Nmin        =   3;
option      =  $1;
mestimator  =  $2;
sigma_s     =  $3;
stdev       =  $4;
EOF

    if (( ${#STRINGS} > 0 )); then
        JEditDetector  -a $DETECTOR  -o $TMPDIR/detector_0.datx  -k "$STRINGS[*]" -q -d 0 >& /dev/null
        JEditDetector  -a $DETECTOR  -o $DETECTOR                -r "$STRINGS[*]" -q -d 0 >& /dev/null
    fi

    PARAMETERS=()

    for ID in ${(k)TRIPODS}; do
        PARAMETERS['JMoveTripod -f tripod.txt -T "$ID addx %" -d 0 >& /dev/null']=$5
        PARAMETERS['JMoveTripod -f tripod.txt -T "$ID addy %" -d 0 >& /dev/null']=$5
        PARAMETERS['JMoveTripod -f tripod.txt -T "$ID addz %" -d 0 >& /dev/null']=$5
    done

    gradient $6 Y

    printf "detector                         %6d %8.4f\n" $N $Y

    if (( ${#STRINGS} > 0 )); then
        JMergeDetector  -a $DETECTOR  -a $TMPDIR/detector_0.datx  -o $DETECTOR >& /dev/null
    fi
}


# -----------------------------------------------------------------------------------------------------
# Stage a.
# This stage can be used to determine the positions of the strings and tripods.
#
# \param  1       fit option
# \param  2       M-estimator
# \param  3       resolution [s]
# \param  4       standard deviations
# \param  5       step size strings [m]
# \param  6       step size tripods [m]
# \param  7       maximum number of extra steps
# -----------------------------------------------------------------------------------------------------
function stage_a()
{
    cat>$TMPDIR/acoustics_fit_parameters.txt<<EOF
`egrep Tmax_s acoustics_fit_parameters.txt`
Nmin        =   3;
option      =  $1;
mestimator  =  $2;
sigma_s     =  $3;
stdev       =  $4;
EOF

    PARAMETERS=()

    for STRING in $STRINGS[*]; do
        PARAMETERS['JMoveDetector -a $DETECTOR -S "$STRING addx %" -d 0 >& /dev/null']=$5
        PARAMETERS['JMoveDetector -a $DETECTOR -S "$STRING addy %" -d 0 >& /dev/null']=$5
        PARAMETERS['JMoveDetector -a $DETECTOR -s "$STRING addz %" -d 0 >& /dev/null']=$5
    done

    for ID in ${(k)TRIPODS}; do
        PARAMETERS['JMoveTripod -f tripod.txt -T "$ID addx %" -d 0 >& /dev/null']=$6
        PARAMETERS['JMoveTripod -f tripod.txt -T "$ID addy %" -d 0 >& /dev/null']=$6
        PARAMETERS['JMoveTripod -f tripod.txt -T "$ID addz %" -d 0 >& /dev/null']=$6
    done

    for STRING in ${HYDROPHONES:|TRANSMITTERS}; do
        PARAMETERS['JEditHydrophone -f hydrophone.txt -S "$STRING rot %" -o hydrophone.txt -q -d 0 >& /dev/null; JEditTransmitter -f transmitter.txt -S "$STRING rot %" -o transmitter.txt -q -d 0 >& /dev/null']=$(($5 / $RADIUS_M))
    done

    for STRING in $TRANSMITTERS[*]; do

        set_variable  MODULE  `getModule -a $DETECTOR -L "$STRING 0"`

        PARAMETERS['JMoveDetector -a $DETECTOR -M "$MODULE addz %" -d 0 >& /dev/null']=$5
        PARAMETERS['JEditHydrophone -f hydrophone.txt -S "$STRING rot %" -o hydrophone.txt -q -d 0 >& /dev/null; JEditTransmitter -f transmitter.txt -S "$STRING rot %" -o transmitter.txt -q -d 0 >& /dev/null']=$(($5 / $RADIUS_M))
    done

    gradient $7 Y

    printf "detector                         %6d %8.4f\n" $N $Y
}


# -----------------------------------------------------------------------------------------------------
# Stage b.
# This stage can be used to determine the stretching and z-positions of individual strings.
#
# \param  1       fit option
# \param  2       M-estimator
# \param  3       resolution [s]
# \param  4       standard deviations
# \param  5       step size stretching [unit]
# \param  6       step size z-position [m]
# \param  7       maximum number of extra steps
# -----------------------------------------------------------------------------------------------------
function stage_b()
{
    cat>$TMPDIR/acoustics_fit_parameters.txt<<EOF
`egrep Tmax_s acoustics_fit_parameters.txt`
Nmin        =   3;
option      =  $1;
mestimator  =  $2;
sigma_s     =  $3;
stdev       =  $4;
EOF

    for STRING in $STRINGS[*]; do

        JEditDetector  -a $DETECTOR  -o $TMPDIR/detector_1.datx  -r "$STRING" -q -d 0 >& /dev/null
        JEditDetector  -a $DETECTOR  -o $DETECTOR                -k "$STRING" -q -d 0 >& /dev/null

        PARAMETERS=()

        PARAMETERS['JMoveDetector -a $DETECTOR -s "$STRING mul %"  -d 0 >& /dev/null; JMoveDetector -a $DETECTOR -s "$STRING subz \$\(\(% * $CAN_ZMAX_M\)\)"  -d 0 >& /dev/null']=$5
        PARAMETERS['JMoveDetector -a $DETECTOR -s "$STRING addz %" -d 0 >& /dev/null']=$6

        gradient $7 Y

        printf "string     %04d                  %6d %8.4f\n" $STRING $N $Y

        JMergeDetector  -a $DETECTOR  -a $TMPDIR/detector_1.datx  -o $DETECTOR >& /dev/null
    done
}


# -----------------------------------------------------------------------------------------------------
# Stage c.
# This stage can be used to determine the z-positions of the modules.
#
# \param  1       fit option
# \param  2       M-estimator
# \param  3       resolution [s]
# \param  4       standard deviations
# \param  5       step size  [m]
# \param  6       maximum number of extra steps
# -----------------------------------------------------------------------------------------------------
function stage_c()
{
    cat>$TMPDIR/acoustics_fit_parameters.txt<<EOF
`egrep Tmax_s acoustics_fit_parameters.txt`
Nmin        =   3;
option      =  $1;
mestimator  =  $2;
sigma_s     =  $3;
stdev       =  $4;
EOF

    PARAMETERS=()

    for STRING in $STRINGS[*]; do       

        for (( FLOOR = 1; $FLOOR <= 18; FLOOR += 1 )); do

            set_variable  MODULE  `getModule -a $DETECTOR -L "$STRING $FLOOR"`

            PARAMETERS['JMoveDetector -a $DETECTOR -M "$MODULE addz %" -d 0 >& /dev/null']=$5
        done
    done
        
    gradient $6 Y
        
    printf "detector                         %6d %8.4f\n" $N $Y
}


# -----------------------------------------------------------------------------------------------------
# Stage d.
# This stage can be used to determine the z-positions of the anchors
#
# \param  1       fit option
# \param  2       M-estimator
# \param  3       resolution [s]
# \param  4       standard deviations
# \param  5       step size  [m]
# \param  6       maximum number of extra steps
# -----------------------------------------------------------------------------------------------------
function stage_d()
{
    cat>$TMPDIR/acoustics_fit_parameters.txt<<EOF
`egrep Tmax_s acoustics_fit_parameters.txt`
Nmin        =   3;
option      =  $1;
mestimator  =  $2;
sigma_s     =  $3;
stdev       =  $4;
EOF

    PARAMETERS=()

    for STRING in $STRINGS[*]; do       

        set_variable  MODULE  `getModule -a $DETECTOR -L "$STRING 0"`

        PARAMETERS['JMoveDetector -a $DETECTOR -M "$MODULE addz %" -d 0 >& /dev/null']=$5
    done
        
    gradient $6 Y
        
    printf "detector                         %6d %8.4f\n" $N $Y
}


# -----------------------------------------------------------------------------------------------------
# Stage D.
# This stage can be used to determine the rotations of the strings and z-positions of the anchors
# including the response of hydrophones.
# It includes a determination of the positions of the strings and tripods.
#
# \param  1       fit option
# \param  2       M-estimator
# \param  3       resolution [s]
# \param  4       standard deviations
# \param  5       step size strings [m]
# \param  6       step size tripods [m]
# \param  7       maximum number of extra steps
# -----------------------------------------------------------------------------------------------------
function stage_D()
{
    cat>$TMPDIR/acoustics_fit_parameters.txt<<EOF
`egrep Tmax_s acoustics_fit_parameters.txt`
Nmin        =   3;
option      =  $1;
mestimator  =  $2;
sigma_s     =  $3;
stdev       =  $4;
EOF

    set_variable  HYDROPHONE                 

    PARAMETERS=()

    for STRING in $STRINGS[*]; do

        set_variable  MODULE  `getModule -a $DETECTOR -L "$STRING 0"`

        PARAMETERS['JMoveDetector -a $DETECTOR -S "$STRING addx %" -d 0 >& /dev/null']=$5
        PARAMETERS['JMoveDetector -a $DETECTOR -S "$STRING addy %" -d 0 >& /dev/null']=$5
        PARAMETERS['JMoveDetector -a $DETECTOR -s "$STRING addz %" -d 0 >& /dev/null']=$5
        PARAMETERS['JMoveDetector -a $DETECTOR -M "$MODULE addz %" -d 0 >& /dev/null']=$5

        PARAMETERS['JEditHydrophone -f hydrophone.txt -S "$STRING rot %" -o hydrophone.txt -q -d 0 >& /dev/null; JEditTransmitter -f transmitter.txt -S "$STRING rot %" -o transmitter.txt -q -d 0 >& /dev/null']=$(($5 / $RADIUS_M))
    done

    for ID in ${(k)TRIPODS}; do
        PARAMETERS['JMoveTripod -f tripod.txt -T "$ID addx %" -d 0 >& /dev/null']=$6
        PARAMETERS['JMoveTripod -f tripod.txt -T "$ID addy %" -d 0 >& /dev/null']=$6
        PARAMETERS['JMoveTripod -f tripod.txt -T "$ID addz %" -d 0 >& /dev/null']=$6
    done

    gradient $7 Y

    printf "detector                         %6d %8.4f\n" $N $Y
}