JEvtReweightGSeaGenParameterScan.sh

Go to the documentation of this file.

#!/bin/zsh
#
# \author bjung
#
version=1.0
script=${0##*/}

# ------------------------------------------------------------------------------------------
#
#                         Utility script to evaluate gSeaGen reweightings
#
# ------------------------------------------------------------------------------------------


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

source $JPP_DIR/setenv.sh $JPP_DIR
source $JPP_DIR/software/scripts/qlib.sh


set_variable     IFS0                     ${IFS}

set_variable:    DEBUG                    REWEIGHTING_DEBUG                           2

set_variable+    BATCH_PROCESSING         REWEIGHTING_BATCH_PROCESSING                1
set_variable+    BATCH_LOGGING            REWEIGHTING_BATCH_LOGGING                   1
set_variable:    QUEUE_SIZE               REWEIGHTING_QUEUE_SIZE                      3
set_variable:    QUEUE_NAME               REWEIGHTING_QUEUE_NAME                      ${TMPDIR:-/tmp}/queue
set_variable:    JOB_SIZE                 REWEIGHTING_JOB_SIZE                        1

set_variable:    TEST_TYPE1D              REWEIGHTING_TEST_TYPE1D                     0
set_variable:    TEST_TYPE1D_PARAMETERS   REWEIGHTING_TEST_TYPE1D_PARAMETERS          0.05
set_variable:    TEST_TYPE2D              REWEIGHTING_TEST_TYPE2D                     2
set_variable:    TEST_TYPE2D_PARAMETERS   REWEIGHTING_TEST_TYPE2D_PARAMETERS          0.05

set_variable:    SCRATCHDIR               REWEIGHTING_SCRATCHDIR                      ${TMPDIR:-/tmp}

set_variable:    N_TRIGGERED_HITS_RANGE   REWEIGHTING_N_TRIGGERED_HITS_RANGE          "20 $((2**31))"
set_variable:    COSTH_RANGE              REWEIGHTING_COSTH_RANGE                     "-1.0 1.0"
set_variable:    ENERGY_RANGE             REWEIGHTING_ENERGY_RANGE                    "0.0 $((2.0**31))"

set_variable:    RECOTYPE                 REWEIGHTING_RECOTYPE                        "track"

set_variable     OSCPARAMETERS                                                        ""
set_variable     WORKDIR                                                              ${TMPDIR:-/tmp}


if do_usage $*; then
    usage "$script <detector file> <input file>+ <particle PDG code <reweighting formula <parameter file>>>+ <oscillation probability table> [oscillation parameter file [working directory]]\n" \
          "Note that if more than one input file is given all optional parameters need to be specified."
fi

case $# in
    8) set_variable WORKDIR                  $8;&
    7) set_variable OSCPARAMETERS            $7;&
    6) set_variable OSCPROBTABLE             $6;
       set_array    PARAMETER_FILES          $5;
       set_array    FORMULAE                 $4;
       set_array    NUTYPES                  $3;
       set_array    INPUT_FILES              $2;
       set_variable DETX_FILE                $1;;
    *) set_variable DETX_FILE                ${argv[1]}
       
       i=2;
       
       while [[ -f ${argv[$i]} ]]; do
           set_array INPUT_FILES \
                     ${INPUT_FILES[@]}       ${argv[$i]}
           (( ++i ))
       done

       while (( i < $# - 2 )); do
           set_array NUTYPES \
                     ${NUTYPES[@]}           ${argv[$i]}
           set_array FORMULAE \
                     ${FORMULAE[@]}          ${argv[$((i+1))]}
           set_array PARAMETER_FILES \
                     ${PARAMETER_FILES[@]}   ${argv[$((i+2))]}
           (( i+=3 ))
       done

       set_variable  OSCPROBTABLE            ${argv[-3]}
       set_variable  OSCPARAMETERS           ${argv[-2]}
       set_variable  WORKDIR                 ${argv[-1]}
esac


# Set neutrino labels

typeset -A NULABELS
set_array  NULABELS \
    -16 "anutau" \
    -14 "anumu"  \
    -12 "anue"   \
     12 "nue"    \
     14 "numu"   \
     16 "nutau"


# Initialise local batch processor

if [[ ${BATCH_PROCESSING} == 1 ]]; then

    if   ( is_CCNikhef ); then

        echo "Submitting reweighting and histogram comparison jobs to nikhef-stbc batch queue: ${QUEUE_NAME}..."

        if [[ ${SCRATCHDIR} == "${TMPDIR:-/tmp}" ]]; then
            fatal "Please specify a valid directory for local \'scratch\' files other than \"${TMPDIR:-/tmp}\""
        fi
        
    elif ( is_CCLyon ); then

        fatal "No batch-processing option available yet for CC-Lyon."
        
    else
        
        echo "Submitting reweighting and histogram comparison jobs to local batch queue: ${QUEUE_NAME}..."
        
        source qlib.sh

        $JPP_DIR/examples/scripts/mkqueue.sh  ${QUEUE_NAME}  ${QUEUE_SIZE}  "start"
    fi
fi


# Retrieve MC files

set_array    MC_MERGE_FILES               $(find ${WORKDIR} -maxdepth 1 -name "*.merge.root")

if (( ${#MC_MERGE_FILES} != 0 )); then

    echo  "----------------------------"    
    printf '%s\n' ${MC_MERGE_FILES[@]}
    echo  "----------------------------"

    echo  "Remove the above regular files? (y/n)"

    read   REMOVE_MC_MERGE_FILES
    
    if [[ ${REMOVE_MC_MERGE_FILES} == "y"   || ${REMOVE_MC_MERGE_FILES} == "Y"   ||
          ${REMOVE_MC_MERGE_FILES} == "yes" || ${REMOVE_MC_MERGE_FILES} == "Yes" || ${REMOVE_MC_MERGE_FILES} == "YES" ]]; then
        
        rm ${MC_MERGE_FILES[@]}
    fi
fi


# Divide into DAQ- and MC-files

echo "Dividing into DAQ- and MC-files..."

for INPUT_FILE in ${INPUT_FILES[@]}; do

    if [[ "$(JPrintHeader -f ${INPUT_FILE} -k simul --)" != " 0" ]]; then # MC-file
        set_array INPUT_MC_FILES          ${INPUT_MC_FILES[@]}  ${INPUT_FILE}
    else                                                                  # DAQ-file
        set_array INPUT_DAQ_FILES         ${INPUT_DAQ_FILES[@]} ${INPUT_FILE}
    fi
done

if (( ${#INPUT_MC_FILES} <= 0 )); then
    fatal "No MC input files specified!"
fi


# Merge MC-files

if [[ ${#MC_MERGE_FILES} == 0 ]]; then
    
    $JPP_DIR/examples/JAAnet/JSortEvt \
        -f "${INPUT_MC_FILES[*]}"                                                          \
        -o  ${WORKDIR}/\%.merge.root                                                       \
        -d  ${DEBUG} --!
fi

set_array MC_MERGE_FILES                  $(find ${WORKDIR} -maxdepth 1 -name "*.merge.root")


# Process reference files

set_variable REFERENCE_POSTFIT_FILE       ${WORKDIR}/reference.postfit.root

if ( ! reuse_file ${REFERENCE_POSTFIT_FILE} ); then

    if (( ${#INPUT_DAQ_FILES} > 0 )); then
        set_array INPUT_REFERENCE_FILES   ${INPUT_DAQ_FILES[@]}
    else
        set_array INPUT_REFERENCE_FILES   ${MC_MERGE_FILES[@]}
    fi

    configure_job \
        ${SCRATCHDIR}/${REFERENCE_POSTFIT_FILE:t:r}.\%.sh 1

    cat >> ${EXEC_FILE} <<EOF
#!/bin/zsh

if [ -z \${JPP_DIR} ]; then echo \"Variable JPP_DIR undefined.\"; exit; fi
 
source \${JPP_DIR}/setenv.sh \${JPP_DIR}

timer_start

JAAPostfit    \\
    -f "${INPUT_REFERENCE_FILES[*]}"                                                   \\
    -o  ${REFERENCE_POSTFIT_FILE}                                                      \\
    -a  ${DETX_FILE}                                                                   \\
    -#energyRange="${ENERGY_RANGE}"                                                    \\
    -#coszenithRange="${COSTH_RANGE}"                                                  \\
    -#triggeredHitsRange="${N_TRIGGERED_HITS_RANGE}"                                   \\
    -R ${RECOTYPE} -W                                                                  \\
    -d  ${DEBUG} --! ${BATCH_LOGGING:+&>> ${LOG_FILE} 2>&1}

timer_stop
timer_print

rm -f ${EXEC_FILE}
EOF

    submit_job \
        ${EXEC_FILE} ${BATCH_PROCESSING:+${QUEUE_NAME} ${SCRATCHDIR}}

    reset_job_config
fi


# Create histogram comparison ASCII steering file

set_array NAMES                           $(JPrintName -f ${REFERENCE_POSTFIT_FILE}:.)

set_variable STEERING_FILE                ${WORKDIR}/steeringfile.txt

rm -f ${STEERING_FILE}

unset HISTOGRAMS1D
unset HISTOGRAMS2D

for NAME in ${NAMES[@]}; do

    set_variable BINFRAC                  $(JPrintResult \
                                                -f ${REFERENCE_POSTFIT_FILE}:^${NAME}\$         \
                                                -F "GetNcells / ((GetNbinsX+2)*(GetNbinsY+2))"  \
                                                -O "%f" --)
    
    if [[ ${NAME} != "META" && ${NAME} != "TObject" ]]; then

        if   (( BINFRAC > 1.0 )); then
            warning "Cannot perform comparison test for histogram ${NAME} with >2 dimensions; skip."
        elif (( BINFRAC < 1.0 )); then
            echo "${NAME}\$   ${TEST_TYPE1D}   ${TEST_TYPE1D_PARAMETERS}" >> ${STEERING_FILE}
            set_array HISTOGRAMS1D        ${HISTOGRAMS1D[@]} ${NAME}
        else
            echo "${NAME}\$   ${TEST_TYPE2D}   ${TEST_TYPE2D_PARAMETERS}" >> ${STEERING_FILE}
            set_array HISTOGRAMS2D        ${HISTOGRAMS2D[@]} ${NAME}        
        fi
    fi
done


# Process MC-files with standard weights

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

    set_variable MC_MERGE_FILE   ${MC_MERGE_FILES[$i]}
    set_variable POSTFIT_FILE    ${MC_MERGE_FILE:r}.postfit.root
    set_variable COMPARISON_FILE ${POSTFIT_FILE:r}.comparison.root

    echo "Submitting histogram comparison job for ${MC_MERGE_FILE:t}"   

    configure_job \
        ${SCRATCHDIR}/${COMPARISON_FILE:t:r}.job\%.sh ${JOB_SIZE}

    cat >> ${EXEC_FILE} <<EOF
${START_JOB+#!/bin/zsh}

${START_JOB+"if [ -z \${JPP_DIR} ]; then echo \"Variable JPP_DIR undefined.\"; exit; fi"}

${START_JOB+"source \${JPP_DIR}/setenv.sh \${JPP_DIR}"}

timer_start

JAAPostfit     \\
    -f ${MC_MERGE_FILE}                                                                    \\
    -o ${POSTFIT_FILE}                                                                     \\
    -a ${DETX_FILE}                                                                        \\
    -@energyRange="${ENERGY_RANGE}"                                                        \\
    -@coszenithRange="${COSTH_RANGE}"                                                      \\
    -@triggeredHitsRange="${N_TRIGGERED_HITS_RANGE}"                                       \\
    -R ${RECOTYPE} -W                                                                      \\
    -d ${DEBUG} --! ${BATCH_LOGGING:+&> ${LOG_FILE}}

timer_stop
timer_print ${BATCH_LOGGING:+&>> ${LOG_FILE} 2>&1}
        
timer_start
        
JZebraMantis       \\
    -a ${REFERENCE_POSTFIT_FILE}                                                           \\
    -b ${POSTFIT_FILE}                                                                     \\
    -s ${STEERING_FILE}                                                                    \\
    -o ${COMPARISON_FILE}                                                                  \\
    -t ${COMPARISON_FILE:r}.txt                                                            \\
    -d ${DEBUG} --! ${BATCH_LOGGING:+&>> ${LOG_FILE} 2>&1}

timer_stop
timer_print ${BATCH_LOGGING:+&>> ${LOG_FILE} 2>&1}

${SUBMIT_JOB+rm -f ${EXEC_FILE}}
EOF

    if (( ${+SUBMIT_JOB} || i == ${#MC_MERGE_FILES} )) ; then
        submit_job \
            ${EXEC_FILE} ${BATCH_PROCESSING:+${QUEUE_NAME} ${SCRATCHDIR} ${JOBIDS[1]}}
    fi
done

reset_job_config


# Process MC-files with redefined weights

set_variable PARDELIM     ';'
set_variable PARSET_LIST  ${WORKDIR}/parametersets.txt

paste -d"$PARDELIM" ${PARAMETER_FILES[@]} > ${PARSET_LIST}

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

    set_variable MC_MERGE_FILE   ${MC_MERGE_FILES[$i]}

    echo "Submitting reweighting job for ${MC_MERGE_FILE:t}"    

    if [[ -f ${MC_MERGE_FILE:r}.reweighting.comparison.meta.txt ]]; then
        rm -i ${MC_MERGE_FILE:r}.reweighting.comparison.meta.txt
    fi

    set_array PRIMARIES \
              $(JPrintHeader -f ${MC_MERGE_FILE} -k primary -- | cut -d'=' -f2)
    
    set_array PRIMARIES \
              ${PRIMARIES[@]} \
              $(JPrintHeader -f ${MC_MERGE_FILE} -k flux | cut -d' ' -f3)
    
    while IFS="$PARDELIM" read -r PARAMETERSETS; do

        unset PARINPUTS

        if [[ ${PARAMETERSETS:0:1} == \# ]]; then continue; fi
        
        for (( j=1; j<=${#PARAMETER_FILES}; ++j )); do

            # Read formula and corresponding parameter set
            
            set_variable NUTYPE           ${NUTYPES[$j]}
            set_variable FORMULA          ${FORMULAE[$j]}
            set_array    PARAMETERS       $(echo $PARAMETERSETS | cut -d"$PARDELIM" -f$j)

            # Check input

            set_variable N_FPAR           $(echo -n $FORMULA | grep -oP "(^|[^x])\[[0-9]*\]" | cut -d'[' -f2- | sort -u | wc -l)
            
            if (( ${#PARAMETERS} != N_FPAR )); then
                fatal "Number of input parameters (= ${N_PARS}) is incompatible with given formula: ${FORMULA}"
            fi

            unset PAREQ

            for (( k=0; k<${#PARAMETERS}; ++k )); do
                set_array PAREQS \
                    ${PAREQS[@]} "p${k}=${PARAMETERS[$((k+1))]};"
            done
            
            for PRIMARY in ${PRIMARIES[@]}; do
                
                if (( $PRIMARY == $NUTYPE )); then

                    set_array PARINPUTS \
                        ${PARINPUTS[@]} "${NULABELS[${NUTYPE}]}=\"${FORMULA}; ${PAREQS[@]}\""
                    break;
                fi
            done
        done

        if (( ${#PARINPUTS} == 0 )); then continue; fi
        
        set_variable REWEIGHTED_FILE      ${MC_MERGE_FILE:r}.reweighting.${i}.root
        set_variable POSTFIT_FILE         ${REWEIGHTED_FILE:r}.postfit.root
        set_variable COMPARISON_FILE      ${POSTFIT_FILE:r}.comparison.root

        configure_job \
            ${SCRATCHDIR}/${COMPARISON_FILE:t:r}.job\%.sh ${JOB_SIZE}
        
        cat >> ${EXEC_FILE} <<EOF
${START_JOB+#!/bin/zsh}

${START_JOB+"if [ -z \${JPP_DIR} ]; then echo \"Variable JPP_DIR undefined.\"; exit; fi"}

${START_JOB+"source \${JPP_DIR}/setenv.sh \${JPP_DIR}"}

timer_start

JEvtReweightGSeaGen \\
    -f  ${MC_MERGE_FILE}                                                                   \\
    -o  ${REWEIGHTED_FILE}                                                                 \\
    -P  ${OSCPROBTABLE}                                                                    \\
    $(printf '-#%s\n' "${PARINPUTS[*]}")                                                   \\
    ${OSCPARAMETERS:+$(printf '-@%s' ${OSCPARAMETERS})}                                    \\
    -d  ${DEBUG} --! ${BATCH_LOGGING:+&> ${LOG_FILE} 2>&1}

timer_stop
timer_print ${BATCH_LOGGING:+&>> ${LOG_FILE} 2>&1}

timer_start

JAAPostfit          \\
    -f ${REWEIGHTED_FILE}                                                                  \\
    -o ${POSTFIT_FILE}                                                                     \\
    -a ${DETX_FILE}                                                                        \\
    -#energyRange="${ENERGY_RANGE}"                                                        \\
    -#coszenithRange="${COSTH_RANGE}"                                                      \\
    -#triggeredHitsRange="${N_TRIGGERED_HITS_RANGE}"                                       \\
    -R ${RECOTYPE} -W                                                                      \\
    -d ${DEBUG} --! ${BATCH_LOGGING:+&>> ${LOG_FILE} 2>&1}

timer_stop
timer_print ${BATCH_LOGGING:+&>> ${LOG_FILE} 2>&1}

timer_start

JZebraMantis        \\
    -a ${REFERENCE_POSTFIT_FILE}                                                           \\
    -b ${POSTFIT_FILE}                                                                     \\
    -s ${STEERING_FILE}                                                                    \\
    -o ${COMPARISON_FILE}                                                                  \\
    -t ${COMPARISON_FILE:r}.txt                                                            \\
    -d ${DEBUG} --! ${BATCH_LOGGING:+&>> ${LOG_FILE} 2>&1}

timer_stop
timer_print ${BATCH_LOGGING:+&>> ${LOG_FILE} 2>&1}

rm -f ${REWEIGHTED_FILE}
${SUBMIT_JOB+rm -f ${EXEC_FILE}}
EOF

        if (( ${+SUBMIT_JOB} || i == ${#MC_MERGE_FILES} )) ; then
            submit_job \
                ${EXEC_FILE} ${BATCH_PROCESSING:+${QUEUE_NAME} ${SCRATCHDIR} ${JOBIDS[1]}}
        fi

        echo "${COMPARISON_FILE:r}.txt; ${PARINPUTS[*]}" >> ${WORKDIR}/comparisonResults.meta.txt
        
    done < ${PARSET_LIST}
done

if (( ${+JOBIDS} )); then
    rm -f ${SCRATCHDIR}/reweighting.jobIDs.txt
    printf '%s\n' ${JOBIDS[@]} > ${SCRATCHDIR}/reweighting.jobIDs.txt
fi

reset_job_config

set_variable IFS                                  ${IFS0}