feat(qa): QA of RG-A Spring 2018 6.4 GeV data

bin/qtl

@@ -18,6 +18,7 @@ usage() {
       physics       Generate and analyze physics QA timelines (Step 2)
       error         Scan for errors in Slurm logs (for Step 1)
       reheat        Reproduce a data file, e.g., to rerun postprocessing
+      xcharge       Cross check and analyze FC charge
       xtrain        Cross check run list from trains and DSTs
 
   OPTIONS: Each command has its own set of options; run a command with no
@@ -42,6 +43,7 @@ case $cmd in
   ph*) exec $TIMELINESRC/bin/qtl-physics "$@" ;;
   er*) exec $TIMELINESRC/bin/qtl-error "$@" ;;
   re*) exec $TIMELINESRC/bin/qtl-reheat "$@" ;;
+  xc*) exec $TIMELINESRC/bin/qtl-xcharge "$@" ;;
   xt*) exec $TIMELINESRC/bin/qtl-xtrain "$@" ;;
   -v|--version)
     echo $(mvn -q help:evaluate -Dexpression=project.version -DforceStdout -f $TIMELINESRC/pom.xml || echo "UNKNOWN")

bin/qtl-reheat

@@ -16,6 +16,7 @@ declare -A REHEAT_METHODS=(
 
 # default options
 dataset=train
+coatjava_version=13.3.0 # the first version with rollover correction, used for RG-A Spring 2018 reheating
 declare -A modes
 for key in submit; do
   modes[$key]=false
@@ -35,12 +36,16 @@ usage() {
   REQUIRED OPTIONS:
      -i [INPUT_DIR]    input directory of HIPO files; all HIPO files will be reheated
      -o [OUTPUT_DIR]   output directory
-     -c [COMMAND]      which reheating method to use, one of:
+     -m [METHOD]       which reheating method to use, one of:
 $(for key in "${!REHEAT_METHODS[@]}"; do printf "%24s %-11s %-s\n" "" "$key" "${REHEAT_METHODS[$key]}"; done)
 
   OPTIONAL OPTIONS:
      -d [DATASET]      unique name for this dataset
                        default: $dataset
+     -c [COATJAVA_VER] coatjava version number, either
+                       - version from coatjava env module
+                       - absolute path to coatjava
+                       default: $coatjava_version
      --submit          submit the slurm jobs, rather than just
                        printing the \`sbatch\` command
   """ >&2
@@ -53,16 +58,17 @@ fi
 # parse options
 inputDir=""
 outputDir=""
-cmd=""
-while getopts "i:o:c:d:h-:" opt; do
+mtd=""
+while getopts "i:o:m:d:c:h-:" opt; do
   case $opt in
     i) inputDir=$OPTARG ;;
     o) outputDir=$OPTARG ;;
-    c)
+    m)
       [ -z "${REHEAT_METHODS[$OPTARG]-}" ] && printError "unknown command '$OPTARG'" && exit 100
-      cmd=$OPTARG
+      mtd=$OPTARG
       ;;
     d) dataset=$OPTARG ;;
+    c) coatjava_version=$OPTARG ;;
     h)
       usage
       exit 101
@@ -80,10 +86,7 @@ done
 # check arguments
 [ -z "$inputDir" ] && printError "missing input directory argument" && exit 100
 [ -z "$outputDir" ] && printError "missing output directory argument" && exit 100
-[ -z "$cmd" ] && printError "missing command argument" && exit 100
-
-# make sure coatjava environment is loaded
-[ -z "${COATJAVA-}" ] && printError "COATJAVA environment variable is not set" && exit 100
+[ -z "$mtd" ] && printError "missing method argument" && exit 100
 
 # get list of input files
 [ ! -d $inputDir ] && printError "input directory '$inputDir' not found" && exit 100
@@ -95,7 +98,7 @@ mkdir -p $outputDir
 outputDir=$(realpath $outputDir)
 
 # start job lists
-jobName=reheat.$cmd.$dataset
+jobName=reheat.$mtd.$dataset
 slurmDir=./slurm
 jobList=$slurmDir/job.$jobName.list
 > $jobList
@@ -122,17 +125,40 @@ echo "INPUT FILE:  $inputFile"
 echo "OUTPUT FILE: $outputFile"
 EOF
 
+  # handle coatjava
+  if [[ "$coatjava_version" == /* ]]; then
+    cat >> $jobScript << EOF
+# set coatjava environment
+export COATJAVA=$coatjava_version
+export CLAS12DIR=\$COATJAVA
+echo "COATJAVA:    \$COATJAVA"
+EOF
+  else
+    cat >> $jobScript << EOF
+# load coatjava module
+source /usr/share/Modules/init/bash
+module purge
+module use /scigroup/cvmfs/hallb/clas12/sw/modulefiles
+module load clas12
+module switch coatjava/$coatjava_version
+echo "COATJAVA:    \$COATJAVA"
+EOF
+  fi
+
   # add main command to the job script
-  case $cmd in
+  case $mtd in
 
     rollover)
       mkdir -p $outputDir/tmp
       tmpFile=$outputDir/tmp/$inputBase.hipo
       cat >> $jobScript << EOF
-echo "TMP FILE:    $tmpFile"
+# remove output files, if they exist
+[ -f "$tmpFile" ] && echo "WARNING: clobbering $tmpFile" >&2 && rm -v $tmpFile
+[ -f "$outputFile" ] && echo "WARNING: clobbering $outputFile" >&2 && rm -v $outputFile
 # fix the clock rollover
-$(which rebuild-scalers) -l INFO -c X -o $tmpFile $inputFile -- -Xmx1536m
-$(which postprocess) -l INFO -q 1 -o $outputFile $tmpFile
+echo "TMP FILE:    $tmpFile"
+\$COATJAVA/bin/rebuild-scalers -l INFO -c X -o $tmpFile $inputFile -- -Xmx1536m
+\$COATJAVA/bin/postprocess -l INFO -q 1 -o $outputFile $tmpFile
 rm -v $tmpFile
 EOF
 

bin/qtl-xcharge

@@ -0,0 +1,87 @@
+#!/usr/bin/env bash
+
+set -euo pipefail
+source $(dirname $0)/../libexec/environ.sh
+
+# constants ############################################################
+# methods: name -> description
+declare -A METHODS=(
+  [charge]="run charge analysis (from Bhawani)"
+  [clock]="plot the clock vs timestamp"
+)
+########################################################################
+
+# usage
+suffix=''
+sep="================================================================"
+usage() {
+  echo """
+  Take a closer look at the FC charge
+  NOTE: chefs should not typically need to run this
+
+  $sep
+  USAGE: qtl xcharge [OPTIONS]
+  $sep
+
+  REQUIRED OPTIONS:
+     -i [INPUT_FILE]   input HIPO file
+     -o [OUTPUT_DIR]   output directory
+     -m [METHOD]       which analysis method, one of:
+$(for key in "${!METHODS[@]}"; do printf "%24s %-11s %-s\n" "" "$key" "${METHODS[$key]}"; done)
+
+  OPTIONAL OPTIONS:
+     -s [SUFFIX]       a suffix to append to output files
+                       default: $suffix
+  """ >&2
+}
+if [ $# -lt 1 ]; then
+  usage
+  exit 101
+fi
+
+# parse options
+inputFile=""
+outputDir=""
+mtd=""
+while getopts "i:o:m:s:h-:" opt; do
+  case $opt in
+    i) inputFile=$OPTARG ;;
+    o) outputDir=$OPTARG ;;
+    m)
+      [ -z "${METHODS[$OPTARG]-}" ] && printError "unknown command '$OPTARG'" && exit 100
+      mtd=$OPTARG
+      ;;
+    s) suffix=$OPTARG ;;
+    h)
+      usage
+      exit 101
+      ;;
+    *) exit 100 ;;
+  esac
+done
+
+# check arguments
+[ -z "$inputFile" ] && printError "missing input directory argument" && exit 100
+[ -z "$outputDir" ] && printError "missing output directory argument" && exit 100
+[ -z "$mtd" ] && printError "missing method argument" && exit 100
+
+# make output directory
+mkdir -p $outputDir
+
+# run the script
+case $mtd in
+  charge)
+    $TIMELINESRC/qa-physics/charge_analysis/analyze_charge.py $inputFile $outputDir $suffix
+    ;;
+  clock)
+    $TIMELINESRC/libexec/run-groovy-timeline.sh $TIMELINESRC/qa-physics/charge_analysis/analyze_clock.groovy $inputFile $outputDir rollover_corr_disabled_$suffix 0
+    $TIMELINESRC/libexec/run-groovy-timeline.sh $TIMELINESRC/qa-physics/charge_analysis/analyze_clock.groovy $inputFile $outputDir rollover_corr_enabled_$suffix  1
+    root -b -q $TIMELINESRC/qa-physics/charge_analysis/plot_clock.C'("'$outputDir'", "'rollover_corr_disabled_$suffix'",0)'
+    root -b -q $TIMELINESRC/qa-physics/charge_analysis/plot_clock.C'("'$outputDir'", "'rollover_corr_enabled_$suffix'",1)'
+    ;;
+  *)
+    echo "ERROR: unknown method '$mtd'" >&2
+    exit 100
+    ;;
+esac
+echo "Done, see files in '$outputDir' with suffix '$suffix'"

qa-physics/charge_analysis/README.md

@@ -1,10 +1,13 @@
+These scripts are called by `qtl xcharge`
+
 # Charge Analysis
 
 From Bhawani Singh
 
-```bash
-./analyze HIPO_FILE
-```
-Produces PNG files comparing the DAQ-gated FC charge determined from
-- directly from the `RUN::scaler` bank
-- from the livetime and ungated charge, by multiplication
+`./analyze_charge.py` Produces PNG files comparing the DAQ-gated FC charge determined from
+- the `RUN::scaler` bank directly
+- the livetime and ungated charge, by multiplication
+
+# Clock analysis
+- `analyze_clock.groovy`: gets the clock for each timestamp
+- `plot_clock.C`: plots the clock vs. timestamp

qa-physics/charge_analysis/analyze.py → qa-physics/charge_analysis/analyze_charge.py

@@ -22,10 +22,18 @@
 
 def main():
 
-    if len(sys.argv) != 2:
-        print(f'USAGE: {sys.argv[0]} [HIPO_FILE]')
+    if len(sys.argv) != 4:
+        print(f'''
+USAGE: {sys.argv[0]} [INPUT_HIPO_FILE] [OUTPUT_DIR] [OUTPUT_FILE_SUFFIX]
+    INPUT_HIPO_FILE       input HIPO file
+    OUTPUT_FILE_SUFFIX    append this string to the output
+                          file name; useful if you are comparing
+                          output files before and after reheating
+        ''')
         exit(2)
-    hipo_file = sys.argv[1]
+    hipo_file     = sys.argv[1]
+    output_dir    = sys.argv[2]
+    output_suffix = sys.argv[3]
 
     hipo_prefix = os.getenv('HIPO')
     if hipo_prefix == None:
@@ -102,17 +110,18 @@ def main():
         ax.tick_params(axis='both', labelsize=9)
 
     fig1.tight_layout(rect=[0, 0.03, 1, 0.95])
-    fig1.savefig(f'fcup_vs_timestamp_{run_number}.png', bbox_inches='tight', dpi=300)
+    fig1.savefig(f'{output_dir}/fcup_vs_timestamp_{run_number}_{output_suffix}.png', bbox_inches='tight', dpi=300)
     plt.close(fig1)
     # ---------- Compute Chunked FCUP Gated with neighbor handling ----------
     chunk_size = 2000
     num_chunks = len(timestamps) // chunk_size
 
-    chunk_caseA, chunk_caseB, chunk_caseC, chunk_default = [], [], [], []
-    cum_caseA, cum_caseB, cum_caseC, cum_default = [], [], [], []
+
+    chunk_caseA, chunk_caseB, chunk_caseC, chunk_default, chunk_default_ungated = [], [], [], [], []
+    cum_caseA, cum_caseB, cum_caseC, cum_default, cum_default_ungated = [], [], [], [], []
     chunk_indices, skipped_counts = [], []
 
-    runA, runB, runC, runDef = 0, 0, 0, 0
+    runA, runB, runC, runDef, runDefUng = 0, 0, 0, 0, 0
     total_skipped = 0
 
     corrected_livetimes_A = []
@@ -130,7 +139,7 @@ def main():
         fcupgated_diff = np.diff(fcupgateds[start:end])
         live_sub = live_times[start+1:end]
 
-        sumA, sumB, sumC, sumDef = 0, 0, 0, 0
+        sumA, sumB, sumC, sumDef, sumDefUng = 0, 0, 0, 0, 0
         skipped_in_chunk = 0
 
         for j, lt in enumerate(live_sub):
@@ -146,6 +155,8 @@ def main():
                 corrected_livetimes_C.append(lt)
                 # Default
                 sumDef += fcupgated_diff[j]
+                # Default ungated
+                sumDefUng += fcup_diff[j]
             else:
                 # ----- Case A/B nearest-neighbor substitution -----
                 idx_candidates = []
@@ -171,6 +182,9 @@ def main():
 
                     # Default
                     sumDef += fcupgated_diff[j]
+
+                    # Default ungated
+                    sumDefUng += fcup_diff[j]
                 else:
                     skipped_in_chunk += 1
                     total_skipped += 1
@@ -192,25 +206,28 @@ def main():
         runB += sumB
         runC += sumC
         runDef += sumDef
+        runDefUng += sumDefUng
 
         chunk_caseA.append(sumA)
         chunk_caseB.append(sumB)
         chunk_caseC.append(sumC)
         chunk_default.append(sumDef)
+        chunk_default_ungated.append(sumDefUng)
         cum_caseA.append(runA)
         cum_caseB.append(runB)
         cum_caseC.append(runC)
         cum_default.append(runDef)
+        cum_default_ungated.append(runDefUng)
         chunk_indices.append(i)
         skipped_counts.append(skipped_in_chunk)
 
     logger.info(f"Computed chunked FCUP Gated values with neighbor handling (Cases A, B, C).")
     logger.info(f"Total skipped events (no valid LT neighbor): {total_skipped}")
 
     # ---------- Plot 2: Chunked FCUP Gated + Ratios + Skips + LT Distribution ----------
-    fig2, (ax_top, ax_mid, ax_bottom, ax_ltdist) = plt.subplots(
-        4, 1, figsize=(12, 14), sharex=False,
-        gridspec_kw={'height_ratios': [3, 1, 1, 2]}
+    fig2, (ax_top, ax_mid, ax_gatedrat, ax_bottom, ax_ltdist) = plt.subplots(
+        5, 1, figsize=(12, 17), sharex=False,
+        gridspec_kw={'height_ratios': [3, 1, 1, 1, 2]}
     )
     fig2.suptitle(f'Run {run_number} - Chunked FCUP Gated (Neighbor Handling)', fontsize=16)
 
@@ -219,6 +236,7 @@ def main():
     #ax_top.plot(chunk_indices, cum_caseB, label='Cumulative Case B (LT_nn × FCUPungated_nn)', color='darkgreen', marker='s')
     ax_top.plot(chunk_indices, cum_caseC, label='Cumulative Case C (20-NN mean × FCUPungated)', color='darkorange', marker='d')
     ax_top.plot(chunk_indices, cum_default, label='Cumulative Default (FCUPgated)', color='blue', marker='^')
+    ax_top.plot(chunk_indices, cum_default_ungated, label='Cumulative Default Ungated (FCUPungated)', color='teal', marker='x',linestyle='--')
     ax_top.set_ylabel('Cumulative Σ', fontsize=11)
     ax_top.grid(True, linestyle='--', alpha=0.6)
     ax_top.legend(fontsize=10)
@@ -228,16 +246,27 @@ def main():
     ratioA = np.divide(cum_caseA, cum_default, out=np.full_like(cum_caseA, np.nan, dtype=float), where=np.array(cum_default) != 0)
     #ratioB = np.divide(cum_caseB, cum_default, out=np.full_like(cum_caseB, np.nan, dtype=float), where=np.array(cum_default) != 0)
     ratioC = np.divide(cum_caseC, cum_default, out=np.full_like(cum_caseC, np.nan, dtype=float), where=np.array(cum_default) != 0)
+    ratioDefUng = np.divide(cum_default_ungated, cum_default, out=np.full_like(cum_default_ungated, np.nan, dtype=float), where=np.array(cum_default) != 0)
 
     ax_mid.plot(chunk_indices, ratioA, label='Case A / Default', color='darkred', marker='o')
     #ax_mid.plot(chunk_indices, ratioB, label='Case B / Default', color='darkgreen', marker='s')
     ax_mid.plot(chunk_indices, ratioC, label='Case C / Default', color='darkorange', marker='d')
+    ax_mid.plot(chunk_indices, ratioDefUng, label='Default FCUP Ungated / Default FCUP gated', color='teal', marker='x',linestyle='--')
     ax_mid.axhline(1.0, color='black', linestyle='--', linewidth=1)
     ax_mid.set_ylabel('Ratio', fontsize=11)
     ax_mid.grid(True, linestyle='--', alpha=0.6)
     ax_mid.legend(fontsize=10)
     ax_mid.tick_params(axis='both', labelsize=10)
 
+    # Gated / Ungated ratio panel
+    ratio_gated_ung = np.divide(cum_default, cum_default_ungated, out=np.full_like(cum_default, np.nan, dtype=float), where=np.array(cum_default_ungated) != 0)
+    ax_gatedrat.plot(chunk_indices, ratio_gated_ung, label='FCUPgated / FCUPungated', color='navy', marker='o')
+    ax_gatedrat.axhline(1.0, color='black', linestyle='--', linewidth=1)
+    ax_gatedrat.set_ylabel('Gated / Ungated', fontsize=11)
+    ax_gatedrat.grid(True, linestyle='--', alpha=0.6)
+    ax_gatedrat.legend(fontsize=10)
+    ax_gatedrat.tick_params(axis='both', labelsize=10)
+
     # Bottom-1: skipped events count
     ax_bottom.bar(chunk_indices, skipped_counts, color='gray', alpha=0.7)
     ax_bottom.set_xlabel(f'Chunk Index (Each = {chunk_size} events)', fontsize=11)
@@ -268,7 +297,7 @@ def main():
     ax_ltdist.tick_params(axis='both', labelsize=10)
 
     fig2.tight_layout(rect=[0, 0.03, 1, 0.95])
-    fig2.savefig(f'chunked_fcupgated_comparison_{run_number}.png', bbox_inches='tight', dpi=300)
+    fig2.savefig(f'{output_dir}/chunked_fcupgated_comparison_{run_number}_{output_suffix}.png', bbox_inches='tight', dpi=300)
     plt.close(fig2)