JEvtReweightMupageParameterScan.sh

Go to the documentation of this file.

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

# ------------------------------------------------------------------------------------------
#
#                         Utility script to evaluate MUPAGE reweightings
#
# ------------------------------------------------------------------------------------------


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

source $JPP_DIR/setenv.sh $JPP_DIR

set_variable     IFS0                     ${IFS}
set_variable     IFS                      $'\n'

set_variable:    DEBUG                    REWEIGHT_MUPAGE_DEBUG                           2

set_variable+    BATCH_PROCESSING         REWEIGHT_MUPAGE_BATCH_PROCESSING                1
set_variable+    BATCH_LOGGING            REWEIGHT_MUPAGE_BATCH_LOGGING                   1
set_variable:    QUEUE_SIZE               REWEIGHT_MUPAGE_QUEUE_SIZE                      3
set_variable:    QUEUE_NAME               REWEIGHT_MUPAGE_QUEUE_NAME                      ${TMPDIR:-/tmp}/queue
set_variable:    JOB_SIZE                 REWEIGHT_MUPAGE_JOB_SIZE                        1

set_variable:    TEST_TYPE1D              REWEIGHT_MUPAGE_TEST_TYPE1D                     0
set_variable:    TEST_TYPE1D_PARAMETERS   REWEIGHT_MUPAGE_TEST_TYPE1D_PARAMETERS          0.05
set_variable:    TEST_TYPE2D              REWEIGHT_MUPAGE_TEST_TYPE2D                     2
set_variable:    TEST_TYPE2D_PARAMETERS   REWEIGHT_MUPAGE_TEST_TYPE2D_PARAMETERS          0.05

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

set_variable:    N_TRIGGERED_HITS_RANGE   REWEIGHT_MUPAGE_N_TRIGGERED_HITS_RANGE          "20 $((2**31))"
set_variable:    COSTH_RANGE              REWEIGHT_MUPAGE_COSTH_RANGE                     "-1.0 1.0"
set_variable:    ENERGY_RANGE             REWEIGHT_MUPAGE_ENERGY_RANGE                    "0.0 $((2.0**31))"

set_variable     WORKDIR                                                                  ${TMPDIR:-/tmp}
set_variable     FORMULA                                                                  "[0]"
set_array        PARAMETERS                                                               "1.0"

if do_usage $*; then
    usage "$script <detector file> <input file>+ [working directory [reweighting formula [parameter setting]+]]"
fi

case $# in
    5) set_variable PARAMETER_FILE        $5;&
    4) set_variable FORMULA               $4;&
    3) set_variable WORKDIR               $3;&
    2) set_array    INPUT_FILES           $2
       set_variable DETX_FILE             $1;;
    *) set_variable DETX_FILE             ${argv[1]}
       set_array    INPUT_FILES           ${argv[2,-4]}
       set_variable WORKDIR               ${argv[-3]}
       set_variable FORMULA               ${argv[-2]}
       set_variable PARAMETER_FILE        ${argv[-1]};;
esac


# Determine number of parameters in given formula

set_variable N_FPAR       $(echo -n ${FORMULA} | grep -oP "(^|[^x])\[[0-9]*\]" | cut -d'[' -f2- | sort -u | wc -l)
set_variable N_PARS       $(( $(wc -w  < ${PARAMETER_FILE}) -    \
                              $(grep "#" ${PARAMETER_FILE} | wc -w) ))
set_variable N_PARSETS    $(( $(wc -l  < ${PARAMETER_FILE}) -    \
                              $(grep "#" ${PARAMETER_FILE} | wc -l) ))

# Check input

if (( N_PARS % N_FPAR != 0 )); then
    fatal "Number of input parameters (= ${N_PARS}) is incompatible with given formula: ${FORMULA}"
fi

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


# 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


# 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

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


# Process DAQ-files

set_variable DAQ_POSTFIT_FILE             ${WORKDIR}/DAQ.postfit.root
    
if ( ! reuse_file ${DAQ_POSTFIT_FILE} ); then

    # Create histograms

    timer_start
    
    JAAMuonPostfit \
        -f "${INPUT_DAQ_FILES[*]}"                                                         \
        -o  ${DAQ_POSTFIT_FILE}                                                            \
        -a  ${DETX_FILE}                                                                   \
        -@energyRange="${ENERGY_RANGE}"                                                    \
        -@coszenithRange="${COSTH_RANGE}"                                                  \
        -@triggeredHitsRange="${N_TRIGGERED_HITS_RANGE}"                                   \
        -W  "standardFlux"                                                                 \
        -d  ${DEBUG} --!

    timer_stop
    timer_print
fi


# Create histogram comparison ASCII steering file

set_array NAMES                           ${(f)$(JPrintName  -f ${DAQ_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 ${DAQ_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


# Merge MC-files

set_variable MC_MERGE_FILE                ${WORKDIR}/\%.merge.root

if [[ ${#MC_MERGE_FILES} == 0 ]]; then
    
    JMergeEvt \
        -f "${INPUT_MC_FILES[*]}"                                                          \
        -o  ${MC_MERGE_FILE}                                                               \
        -d  ${DEBUG} --!
fi

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


# Process MC-files with standard weights

unset STANDARD_OUTPUT

for MC_MERGE_FILE in ${MC_MERGE_FILES[@]}; do

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

    set_variable EXEC_FILE                ${SCRATCHDIR}/${MC_MERGE_FILE:t:r}.sh
    set_variable LOG_FILE                 ${EXEC_FILE:r}.log

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

timer_start

JAAMuonPostfit     \\
    -f ${MC_MERGE_FILE}                                                                    \\
    -o ${MC_MERGE_FILE:r}.postfit.root                                                     \\
    -a ${DETX_FILE}                                                                        \\
    -@energyRange="${ENERGY_RANGE}"                                                        \\
    -@coszenithRange="${COSTH_RANGE}"                                                      \\
    -@triggeredHitsRange="${N_TRIGGERED_HITS_RANGE}"                                       \\
    -W "standardFlux"                                                                      \\
    -d ${DEBUG} --! ${BATCH_LOGGING:+&> ${LOG_FILE}}

timer_stop
timer_print ${BATCH_LOGGING:+&>> ${LOG_FILE} 2>&1}
        
timer_start
        
JZebraMantis       \\
    -a ${DAQ_POSTFIT_FILE}                                                                 \\
    -b ${MC_MERGE_FILE:r}.postfit.root                                                     \\
    -s ${STEERING_FILE}                                                                    \\
    -o ${MC_MERGE_FILE:r}.postfit.comparison.root                                          \\
    -t ${MC_MERGE_FILE:r}.postfit.comparison.txt                                           \\
    -d ${DEBUG} --! ${BATCH_LOGGING:+&>> ${LOG_FILE} 2>&1}

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

rm -f ${EXEC_FILE}
EOF

    chmod +x    ${EXEC_FILE}

    if [[ ${BATCH_PROCESSING} == 1 ]]; then
        
        if ( is_CCNikhef ); then
            
            set_array JOBIDS \
                ${JOBIDS[@]}                                                               \
            ${$(qsub         \
                -N ${EXEC_FILE:t:r}.${USER}                                                \
                -q ${QUEUE_NAME}                                                           \
                -o ${SCRATCHDIR}                                                           \
                -e ${SCRATCHDIR}                                                           \
                -j oe                                                                      \
                ${EXEC_FILE})%%.*}
        else
            
            puts_queue  ${QUEUE_NAME}     ${EXEC_FILE}
        fi
        
    else

        zsh ${EXEC_FILE}
    fi

    set_array STANDARD_OUTPUT             ${MC_MERGE_FILE:r}.postfit.comparison.root
done


# Process MC-files with redefined weights

unset      REWEIGHTED_OUTPUT
typeset -A REWEIGHTED_OUTPUT

for MC_MERGE_FILE in ${MC_MERGE_FILES[@]}; do

    unset   i
    integer i=1

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

    while read PARAMETERS; do

        if [[ ${${(z)PARAMETERS}[1]} == "#" ]]; then continue; fi # skip comments

        if (( ${#${(z)PARAMETERS}} != N_FPAR )); then
            fatal "Parameterset (${PARAMETERS}) is incompatible with given formula (${FORMULA})"
        fi

        unset PARSET

        for (( j=0; j < N_FPAR; ++j )); do
            PARSET+=" p${j}=$(printf '%.2e' ${${(z)PARAMETERS}[$((j+1))]});"
        done
        
        set_variable OUTPUT_FILE          ${MC_MERGE_FILE:r}.reweighting.${i}.root

        if (( i % JOB_SIZE == 1 )); then

            if (( i + JOB_SIZE - 1 <= N_PARSETS )); then
                set_variable EXEC_FILE    ${SCRATCHDIR}/${MC_MERGE_FILE:t:r}.reweighting.${i}.$((i+JOB_SIZE-1)).sh
            else
                set_variable EXEC_FILE    ${SCRATCHDIR}/${MC_MERGE_FILE:t:r}.reweighting.${i}.${N_PARSETS}.sh
            fi
            
            set_variable START_JOB  1
            unset        SUBMIT_JOB

            rm -f ${EXEC_FILE}
            
        elif (( i % JOB_SIZE == 0 || i == N_PARSETS )); then

            set_array JOBRANGE            ${${(s/./)EXEC_FILE}[-3,-2]}
            
            echo "Submitting reweighting jobs ${JOBRANGE[1]} to ${JOBRANGE[2]} for ${MC_MERGE_FILE:t}"
            
            unset        START_JOB
            set_variable SUBMIT_JOB 1
            
        else
            
            unset START_JOB
            unset SUBMIT_JOB
        fi

        set_variable LOG_FILE             ${EXEC_FILE:r}.log
        
        eval 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

JEvtReweightMupage \\
    -f  ${MC_MERGE_FILE}                                                                   \\
    -o  ${OUTPUT_FILE}                                                                     \\
    -F "${FORMULA}"                                                                        \\
    -@ "${PARSET}"                                                                         \\
    -d  ${DEBUG} --! ${BATCH_LOGGING:+&> ${LOG_FILE} 2>&1}

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

JAAMuonPostfit     \\
    -f ${OUTPUT_FILE}                                                                      \\
    -o ${OUTPUT_FILE:r}.postfit.root                                                       \\
    -a ${DETX_FILE}                                                                        \\
    -@energyRange="${ENERGY_RANGE}"                                                        \\
    -@coszenithRange="${COSTH_RANGE}"                                                      \\
    -@triggeredHitsRange="${N_TRIGGERED_HITS_RANGE}"                                       \\
    -W "rescaledFlux"                                                                      \\
    -d ${DEBUG} --! ${BATCH_LOGGING:+&>> ${LOG_FILE} 2>&1}

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

timer_start

JZebraMantis       \\
    -a ${DAQ_POSTFIT_FILE}                                                                 \\
    -b ${OUTPUT_FILE:r}.postfit.root                                                       \\
    -s ${STEERING_FILE}                                                                    \\
    -o ${OUTPUT_FILE:r}.postfit.comparison.root                                            \\
    -t ${OUTPUT_FILE:r}.postfit.comparison.txt                                             \\
    -d ${DEBUG} --! ${BATCH_LOGGING:+&>> ${LOG_FILE} 2>&1}

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

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

        if (( ${+SUBMIT_JOB} )); then

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

                chmod +x ${EXEC_FILE}

                if ( is_CCNikhef ); then
                    
                    set_array JOBIDS \
                        ${JOBIDS[@]}                                                       \
                    ${$(qsub \
                        -N ${EXEC_FILE:t:r}.${USER}                                        \
                        -q ${QUEUE_NAME}                                                   \
                        -o ${SCRATCHDIR}                                                   \
                        -e ${SCRATCHDIR}                                                   \
                        -j oe                                                              \
                        ${EXEC_FILE})%%.*}
                else
                    
                    puts_queue  ${QUEUE_NAME}     ${EXEC_FILE}
                fi
                
            else
                
                zsh ${EXEC_FILE}
            fi
        fi

        echo "${OUTPUT_FILE:r}.postfit.comparison.txt; ${FORMULA}; ${PARSET}" >> ${WORKDIR}/comparisonResults.meta.txt

        i+=1
        
    done < ${PARAMETER_FILE}
done

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

set_variable IFS                                  ${IFS0}