transmittance.sh

Go to the documentation of this file.

#!/usr/bin/env zsh
#
#
# \author mdejong
#
version=1.0
script=${0##*/}
 
# ------------------------------------------------------------------------------------------
#
#                         Utility script to process OMGsim output.
#
# ------------------------------------------------------------------------------------------
 
if [ -z $JPP_DIR ]; then
    echo "Variable JPP_DIR undefined."
    exit
fi
 
source $JPP_DIR/setenv.sh $JPP_DIR
 
zmodload zsh/mathfunc
 
set_variable:  DEBUG       OMGSIM_WORKDIR         2
set_variable:  FACTOR      OMGSIM_FACTOR          0.5
set_variable:  RECYCLING   OMGSIM_RECYCLING       64
 
if do_usage $*; then
    usage "$script <data directory>"
fi
 
if (( $# != 1 )); then
    fatal "Wrong number of arguments."
fi
 
set_variable  DATA                $1
 
set_variable  DETECTOR            $DATA/singleDOM_OMGsim.detx
set_variable  FILE_NAME_PATTERN   $DATA/DOM%_QE2.root
set_variable  OMGSIM              $JPP_DIR/examples/JCalibrate/OMGsim.txt
 
 
function getv() { JPrintProperties -f $OMGSIM -k $1 -@ "cc=#" }
 
typeset -A BUFFER
 
let "BUFFER[K40Glass]   = `getv glass.K40  ` * `getv mass.glass`"
let "BUFFER[Th232Glass] = `getv glass.Th232` * `getv mass.glass`"
let "BUFFER[U238Glass]  = `getv glass.U238 ` * `getv mass.glass`"
let "BUFFER[K40Water]   = `getv water.K40  ` * (4.0/3.0) * acos(-1.0) * (`getv radius.K40 `**3 - `getv radius.glass`**3)"
let "BUFFER[U238Water]  = `getv water.U238 ` * (4.0/3.0) * acos(-1.0) * (`getv radius.U238`**3 - `getv radius.glass`**3)"
 
echo Radio activities:
 
for KEY VALUE in ${(kv)BUFFER}; do
    printf "%-12s %12.3e [Bq]\n" $KEY $VALUE 
done
 
Total=(${(k)BUFFER})
Water=(`echo ${(k)BUFFER} | tr ' ' '\n' | sed -n 's/\‍(.*Water\‍)/\1/p'`)
Glass=(`echo ${(k)BUFFER} | tr ' ' '\n' | sed -n 's/\‍(.*Glass\‍)/\1/p'`)
 
 
cat>PMT_parameters.txt<<EOF
%.TTS_ns           = -1.0
%.gain             =  1.0
%.gainSpread       =  0.4
%.PunderAmplified  =  0.05
%.slewing          =  1
EOF
 
OPTIONS=("BELL_SHAPE = 2.0;")
 
 
for KEY VALUE in ${(kv)BUFFER}; do
 
    set_variable  INPUT_FILE   ${FILE_NAME_PATTERN/\%/$KEY}
    set_variable  OUTPUT_FILE  \[${KEY}\].root
 
    echo "Processing $INPUT_FILE -> $KEY.log"
 
    nohup \
            $JPP_DIR/examples/JCalibrate/JOMGsim  \
            -a $DETECTOR                          \
            -f $INPUT_FILE                        \
            -o $OUTPUT_FILE                       \
            -R $VALUE                             \
            -P PMT_parameters.txt                 \
            -P "QE = $FACTOR"                     \
        -N $RECYCLING                         \
            -d $DEBUG --!                         >&  $KEY.log &
done
 
wait
 
H2=1.2R
 
for Source in Total Water Glass; do
 
    SOURCE=${(U)Source}
 
    JOpera2D \
        -f\[${^${(P)Source}}\].root:$H2           \
        -u Add                                    \
        -O =                                      \
        -o \[$Source\].root                       >& /dev/null
 
    JProject2D \
        -f \[$Source\].root:$H2                   \
        -o \[$Source\].px.root                    \
        -P x                                      \
        -y "-1.0e10 +1.0e10"                      >& /dev/null
 
    JFitK40 \
        -a $DETECTOR                              \
        -f \[$Source\].root                       \
        -@${^OPTIONS}                             \
        -o fit\[$Source\].root                    \
        -w                                        \
        -M                                         > $Source.log
 
    # fitted function
 
    eval $SOURCE=\‍(`sed -n 's/parameters\..*\.set(\‍([^)]*\‍));/\1/p' $Source.log`\‍)
 
done
 
 
for Source in Water Glass; do
 
    SOURCE=${(U)Source}
 
    # ratio of fitted functions
 
    eval ${SOURCE}\[1\]=$((${(P)${SOURCE}[1]} / $TOTAL[1]))
    
    for (( i = 2; $i <= 5; i += 1 )); do
        eval ${SOURCE}\[$i\]=$((${(P)${SOURCE}[$i]} - $TOTAL[$i]))
    done
 
    # histogram
 
    JOpera1D \
        -f \[$Source\].px.root:\.\*    \
        -f \[Total\].px.root:\.\*      \
        -u Divide                      \
        -o R\[$Source\].root           >& /dev/null
 
    # model
 
    $JPP_DIR/examples/JCalibrate/JPlotK40 \
        -@ "R  = ${(P)${SOURCE}[1]}"   \
        -@ "p1 = ${(P)${SOURCE}[2]}"   \
        -@ "p2 = ${(P)${SOURCE}[3]}"   \
        -@ "p3 = ${(P)${SOURCE}[4]}"   \
        -@ "p4 = ${(P)${SOURCE}[5]}"   \
        -o M\[$Source\].root
done
 
 
function typeout()
{
    PARAMETERS=(R p1 p2 p3 p4)
 
    SOURCE=${(U)1}
 
    printf "%s:\n" $SOURCE
 
    for (( i = 1; $i <= ${#PARAMETERS}; i += 1 )); do
        printf "parameters.%-2s.set(%9.6f);\n" $PARAMETERS[$i]  ${(P)${SOURCE}[$i]}
    done
}
 
typeout  Total
typeout  Water
typeout  Glass
 
 
JPlot1D \
    -w 1200x600                        \
    -f R\[Water\].root:\.\*            \
    -f M\[Water\].root:\.\*            \
    -f R\[Glass\].root:\.\*            \
    -f M\[Glass\].root:\.\*            \
    -g 2                               \
    -> "PMT pair"                     \
    -\^ "Ratio"                        \
    -O HIST                            \
    -T ""                              \
    -L TR                              \
    -o transmittance.gif