JARCAEnergyCorrection.sh

Go to the documentation of this file.

#!/usr/bin/env zsh
#
# \author mdejong
#
version=1.0
script=${0##*/}
 
# ------------------------------------------------------------------------------------------
#
#    Utility script to determine muon energy correction base on Monte Carlo simulations.
#
# ------------------------------------------------------------------------------------------
 
if [ -z $JPP_DIR ]; then
    echo "Variable JPP_DIR undefined."
    exit
fi
 
source $JPP_DIR/setenv.sh $JPP_DIR
source $JPP_DIR/externals/km3net-dataformat/definitions/reconstruction.sh
 
zmodload zsh/mathfunc
 
set_variable   DEBUG              2
set_variable   WORKDIR            ./
set_variable:  NUMBER_OF_EVENTS   RECONSTRUCTION_NUMBER_OF_EVENTS   100000
set_variable:  NUMBER_OF_BINS     RECONSTRUCTION_NUMBER_OF_BINS     30
set_variable:  LOGE               RECONSTRUCTION_LOGE               2.5 9.0
set_variable:  FORMAT             GRAPHICS_FORMAT        gif
set_variable+  BATCH              GRAPHICS_BATCH         -B
 
if do_usage $*; then
    usage "$script <detector> <option> [working directory]"\
          "\nPossible options: run fit interpolate"
fi
 
case $# in
    3) set_variable WORKDIR     $3;&
    2) set_variable OPTION      $2;
       set_variable DETECTOR    $1;;
    *) fatal "Wrong number of arguments."
esac
 
if [[ $OPTION != "run" && $OPTION != "fit" && $OPTION != "interpolate" ]]; then
    fatal "Invalid option $OPTION."
fi
 
if [[ ! -d $WORKDIR ]]; then
    mkdir -p  $WORKDIR
fi
 
let "XMIN = ${LOGE%% *}"
let "XMAX = ${LOGE##* }"
 
if [[ $OPTION == "run" ]]; then
 
    if [[ ! -f $WORKDIR/generator.root ]]; then
 
        eval `JPrintDetector -a $DETECTOR -O CAN`
 
        let "RADIUS_M = $CAN_RADIUS_M + 200.0"
        let "ZMIN_M   = 0.0"
        let "ZMAX_M   = $CAN_ZMAX_M   + 200.0"
 
        JGenerator \
            -C "$CAN_X_M $CAN_Y_M $RADIUS_M $ZMIN_M $ZMAX_M"   \
            -E "$LOGE"                                         \
            -P 13                                              \
            -n $NUMBER_OF_EVENTS                               \
            -o $WORKDIR/generator.root                         \
            -d $DEBUG
    fi
 
    if [[ ! -f $WORKDIR/sirene.root ]]; then
 
        JSirene.sh                  $DETECTOR             $WORKDIR/generator.root           $WORKDIR/sirene.root
 
        rm -f  $WORKDIR/generator.root
    fi
 
    if [[ ! -f $WORKDIR/trigger_efficiency.root ]]; then
 
        set_variable  TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY
 
        cat>$WORKDIR/trigger_parameters.txt<<EOF
trigger3DMuon.enabled           =  1;
trigger3DMuon.numberOfHits      =  5;
trigger3DMuon.numberOfModules   =  5;
trigger3DMuon.gridAngle_deg     = 10;
ctMin           =  0.0;
TMaxLocal_ns    = 15.0;
writeSummary    = 100;
EOF
    
        JTriggerEfficiency.sh       $DETECTOR  $DETECTOR  $WORKDIR/sirene.root              $WORKDIR/trigger_efficiency.root  $WORKDIR/trigger_parameters.txt  $JPP_DATA/PMT_parameters.txt
 
        rm -f  $WORKDIR/sirene.root
    fi
 
    if (( 0 )) && [[ ! -f $WORKDIR/jmuon.root ]]; then
 
        set_variable  RECONSTRUCTION_NUMBER_OF_THREADS    10
 
        JARCAMuonReconstruction.sh  $DETECTOR             $WORKDIR/trigger_efficiency.root  $WORKDIR/%.root                   $JPP_DATA/PMT_parameters.txt
 
        rm -f  $WORKDIR/trigger_efficiency.root
        rm -f  $WORKDIR/jmuon.offline.root
    fi
 
    if (( 1 )) && [[ ! -f $WORKDIR/jmuon.root ]]; then
 
        set_variable  PDF                  $JPP_DATA/J%p.dat
        set_variable  NUMBER_OF_THREADS    10
        set_variable  RATE_HZ              6.0e3
 
        typeset -A PARAMETERS
        typeset -a APPLICATIONS
 
        APPLICATIONS=()
 
        JMuonMCEvt \
            -f $WORKDIR/trigger_efficiency.root          \
            -o $WORKDIR/jmcevt.root                      \
            -d $DEBUG --!
 
        rm -f  $WORKDIR/trigger_efficiency.root
 
        APPLICATIONS+=$JMCEVT
 
        PARAMETERS=(
            roadWidth_m      300.0
            R_Hz             $RATE_HZ
            EMin_log         0.0
            EMax_log         8.0
            TMin_ns           -50.0
            TMax_ns          +450.0
            mestimator       1
            numberOfPrefits  1)
 
        eval JMuonEnergy \
             -a $DETECTOR                   \
             -f $WORKDIR/jmcevt.root        \
             -o $WORKDIR/jenergy.root       \
             -P $PDF                        \
             -@ `make_equation PARAMETERS`  \
             -N ${NUMBER_OF_THREADS}        \
             -d $DEBUG --!
 
        rm -f   $WORKDIR/jmcevt.root
 
        APPLICATIONS+=$JMUONENERGY
 
        JHistory \
            -f $WORKDIR/jenergy.root       \
            -o $WORKDIR/jmuon.root         \
            -A "$APPLICATIONS[*]"          \
            -d $DEBUG --!
 
        rm -f  $WORKDIR/jenergy.root
    fi
fi
 
if [[ $OPTION == "fit" || $OPTION == "interpolate" ]]; then
 
    if ( ! reuse_file $WORKDIR/energy_correction.root ); then
        set_variable  INPUT_FILE  $WORKDIR/jmuon.root
    else
        set_variable  INPUT_FILE  $WORKDIR/energy_correction.root
    fi
 
    if   [[ $OPTION == "fit" ]]; then
 
        let "X0 = 0.5 * ($XMIN + $XMAX)"                 # central value
        let "XR = ($XMAX - $XMIN)"                       # range
 
        X="(x - $X0) / $XR"                              # normalisation
 
        P0="(1.0)"                                       # Legendre polynomials
        P1="($X)"
        P2="(0.5*(3*$X*$X - 1.0))"
        P3="(0.5*$X*(5*$X*$X* - 3.0))"
 
        JEnergyCorrection \
            -f $INPUT_FILE                               \
            -o $WORKDIR/energy_correction.root           \
            -F "[0]*$P0 + [1]*$P1 + [2]*$P2 + [3]*$P3"   \
            -@ "p0 = 0.0"                                \
            -@ "p1 = 1.0"                                \
            -@ "p2 = 0.0"                                \
            -@ "p3 = 0.0"                                \
            -N $NUMBER_OF_BINS                           \
            -x "$LOGE"                                   \
            -d $DEBUG --!
 
    elif [[ $OPTION == "interpolate" ]]; then
 
        JEnergyCorrection \
            -f $INPUT_FILE                               \
            -o $WORKDIR/energy_correction.root           \
            -N $NUMBER_OF_BINS                           \
            -x "$LOGE"                                   \
            -d $DEBUG --!
    fi
 
    JLine \
        -@ "color = 1; style = 2"                        \
        -p "$XMIN $XMIN $XMAX $XMAX"                     \
        -o $WORKDIR/line.root
 
    JPlot2D \
        -f  $WORKDIR/energy_correction.root:h2           \
        -f  $WORKDIR/line.root:\.\*                      \
        -> "E_{fit} [GeV]"                              \
        -< "E_{true} [GeV]"                             \
        -XX                                              \
        -YY                                              \
        -Z                                               \
        -O COLZ                                          \
        -T ""                                            \
        -o ExE.$FORMAT                                   $BATCH
 
    JPlot1D \
        -w 1200x600                                      \
        -f $WORKDIR/energy_correction.root:"py\[.*\]"    \
        -XX                                              \
        -> "E [GeV]"                                    \
        -\^ "number of events [a.u]"                     \
        -T ""                                            \
        -o E1.$FORMAT                                    $BATCH
 
    JPlot1D \
        -f  $WORKDIR/energy_correction.root:h1           \
        -f  $WORKDIR/line.root:\.\*                      \
        -y "$LOGE"                                       \
        -XX                                              \
        -YY                                              \
        -> "log(E_{fit}/GeV)"                           \
        -\^ "log(E_{true}/GeV)"                          \
        -G XY                                            \
        -T ""                                            \
        -o E2.$FORMAT                                    $BATCH
 
    rm -f  $WORKDIR/line.root
 
    echo result stored in $WORKDIR/energy_correction.root
fi