Update of scripts

	modified:   .gitignore
	modified:   README.md
	new file:   bin/hc_coin_replay
	new file:   bin/hc_shms_replay
	new file:   bin/latest_run
	modified:   db2/auto_standard.kinematics
	new file:   db2/auto_standard_hms.kinematics
	new file:   db2/auto_standard_shms.kinematics
	new file:   online_bin/start_auto_replay
	new file:   scripts/good_shms_counter.cxx
	modified:   scripts/make_human_table.cxx
	new file:   scripts/make_jpsi_table.cxx
	new file:   scripts/muon_pion_separation.cxx
	new file:   scripts/plot_shms_singles.cxx
	renamed:    scripts/replay_production_coin_sidis.cxx -> scripts/replay_production_coin.cxx
	new file:   scripts/replay_production_shms.cxx

.gitignore

 */*.nfs*
 raw
 raw*
-bin
 DBASE
 analysis
 ROOTfiles

README.md

@@ -5,6 +5,22 @@ CSV online files
 * [`scripts` directory](scripts/README.md)
 * [`online_monitor` directory](online_monitor/README.md)
 
+## COIN/SHMS/HMS replay
+
+```
+./bin/hc_shms_replay -r run_number -n Nevents
+./bin/hc_shms_replay -h
+./bin/hc_coin_replay -r run_number -n Nevents
+```
+
+### Tips
+
+#### Running 10 replays in parallel using gnu parallel
+```
+parallel -N1 -j10 ./bin/hc_shms_replay -r {}  ::: $(seq 7100 7142)
+```
+
+
 ## Quick Start
 
 [![asciicast](https://asciinema.org/a/222886.svg)](https://asciinema.org/a/222886)

bin/hc_coin_replay

+#!/bin/bash
+#set -u 
+
+function print_the_help {
+  echo "USAGE: $0 -r <run_number> -n <num> "
+  echo "  OPTIONS: "
+  echo "            -r,--run           Required run number"
+  echo "            -n,--n-events        Required number of eventsrun number"
+  #echo "            -o,--online-only   Log online only"
+  #echo "            -a,--all           Log all plots (online, detectors, physics)"
+  #echo "            -n,--no-log        Do not submit any log entries"
+  #echo "            -y,--yes           Automatically submit log entries (does not prompt for y/n)"
+  exit 
+}
+
+function yes_or_no {
+  while true; do
+    read -p "$* [y/n]: " yn
+    case $yn in
+      [Yy]*) return 0 ;;
+      [Nn]*) echo "No entered" ; return 1 ;;
+    esac
+  done
+}
+
+run_number=
+n_events=-1
+container=
+
+if [[ $# -eq 0 ]] ; then
+  print_the_help
+  exit 
+fi
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+  key="$1"
+
+  case $key in
+    -h|--help)
+      shift # past argument
+      print_the_help
+      ;;
+    -r|--run)
+      run_number="$2"
+      shift # past argument
+      shift # past value
+      ;;
+    -n|--n-events)
+      n_events="$2"
+      shift # past argument
+      shift # past value
+      ;;
+    #-c|--container)
+    #  container=1
+    #  #shift # past argument
+    #  shift # past value
+    #  ;;
+    #-a|--all)
+    #  LOG_ALL_PLOTS=1
+    #  shift # past argument
+    #  #shift # past value
+    #  ;;
+    #-o|--online-only)
+    #  ONLINE_ONLY=1
+    #  shift # past argument
+    #  #shift # past value
+    #  ;;
+    *)    # unknown option
+      #POSITIONAL+=("$1") # save it in an array for later
+      echo "unknown option"
+      print_the_help
+      shift # past argument
+      ;;
+  esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
+
+if [[ -z ${n_events} ]] ; then
+  echo  " need argument -n <num> "
+fi
+if [[ -z ${run_number} ]] ; then
+  echo  " need argument -r <run> "
+fi
+
+#if [[ $# -eq 0 ]] ; then
+#  print_the_help
+#  exit 
+#fi
+
+source /group/c-csv/local/setup.sh
+
+# For some reason the module files need to put in a specific order.
+# This is probably because we are using the old TCL version.
+module unload csv/dev
+module load csv/dev
+module load tmux/latest
+module load python/latest
+module load python/2.7.15
+module load python3/latest
+module load vim/latest
+module load git/latest
+module load ruby/2.5.3
+
+
+hcana -b -q "scripts/replay_production.cxx(${run_number},${n_events})"
+
+

bin/hc_shms_replay

+#!/bin/bash
+#set -u 
+
+function print_the_help {
+  echo "USAGE: $0 -r <run_number> -n <num> "
+  echo "  OPTIONS: "
+  echo "            -r,--run           Required run number"
+  echo "            -n,--n-events        Required number of eventsrun number"
+  #echo "            -o,--online-only   Log online only"
+  #echo "            -a,--all           Log all plots (online, detectors, physics)"
+  #echo "            -n,--no-log        Do not submit any log entries"
+  #echo "            -y,--yes           Automatically submit log entries (does not prompt for y/n)"
+  exit 
+}
+
+function yes_or_no {
+  while true; do
+    read -p "$* [y/n]: " yn
+    case $yn in
+      [Yy]*) return 0 ;;
+      [Nn]*) echo "No entered" ; return 1 ;;
+    esac
+  done
+}
+
+run_number=
+n_events=-1
+container=
+
+if [[ $# -eq 0 ]] ; then
+  print_the_help
+  exit 
+fi
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+  key="$1"
+
+  case $key in
+    -h|--help)
+      shift # past argument
+      print_the_help
+      ;;
+    -r|--run)
+      run_number="$2"
+      shift # past argument
+      shift # past value
+      ;;
+    -n|--n-events)
+      n_events="$2"
+      shift # past argument
+      shift # past value
+      ;;
+    #-c|--container)
+    #  container=1
+    #  #shift # past argument
+    #  shift # past value
+    #  ;;
+    #-a|--all)
+    #  LOG_ALL_PLOTS=1
+    #  shift # past argument
+    #  #shift # past value
+    #  ;;
+    #-o|--online-only)
+    #  ONLINE_ONLY=1
+    #  shift # past argument
+    #  #shift # past value
+    #  ;;
+    *)    # unknown option
+      #POSITIONAL+=("$1") # save it in an array for later
+      echo "unknown option"
+      print_the_help
+      shift # past argument
+      ;;
+  esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
+
+if [[ -z ${n_events} ]] ; then
+  echo  " need argument -n <num> "
+fi
+if [[ -z ${run_number} ]] ; then
+  echo  " need argument -r <run> "
+fi
+
+#if [[ $# -eq 0 ]] ; then
+#  print_the_help
+#  exit 
+#fi
+
+source /group/c-csv/local/setup.sh
+
+# For some reason the module files need to put in a specific order.
+# This is probably because we are using the old TCL version.
+module unload csv/dev
+module load csv/dev
+module load tmux/latest
+module load python/latest
+module load python/2.7.15
+module load python3/latest
+module load vim/latest
+module load git/latest
+module load ruby/2.5.3
+
+
+hcana -b -q "scripts/replay_production_shms.cxx(${run_number},${n_events})"
+
+

bin/latest_run

+#!/bin/bash
+
+function print_the_help {
+  echo "USAGE: latest_run [--type=coin] "
+  echo "  OPTIONS: "
+  echo "            -t,--type          Run type (coin, hms, shms) [default: coin]"
+  echo "            -d,--dir           raw dir to use [default: /net/cdaq/cdaql3data/coda/data/raw]"
+  echo "            -c,--cache         Use the mss cache directory (same as using --dir /cache/hallc/spring17/raw)"
+  echo "            -h,--help          print help"
+  exit 
+}
+
+function yes_or_no {
+  while true; do
+    read -p "$* [y/n]: " yn
+    case $yn in
+      [Yy]*) return 0 ;;
+      [Nn]*) echo "No entered" ; return 1 ;;
+    esac
+  done
+}
+#yes_or_no "Upload these ? " && some_command
+
+#if [[ $# -eq 0 ]] ; then
+#  print_the_help
+#  exit 
+#fi
+
+CODATYPE=coin
+USE_RAW_DIR=
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+  key="$1"
+
+  case $key in
+    -h|--help)
+      shift # past argument
+      print_the_help
+      ;;
+    -t|--type)
+      CODATYPE="$2"
+      shift # past argument
+      shift # past value
+      ;;
+    -d|--dir)
+      USE_RAW_DIR="$2"
+      
+      shift # past argument
+      shift # past value
+      ;;
+    -c|--cache)
+      USE_RAW_DIR="/cache/hallc/spring17/raw"
+      shift # past value
+      ;;
+    #-a|--all)
+    #  LOG_ALL_PLOTS=1
+    #  shift # past argument
+    #  #shift # past value
+    #  ;;
+    #-o|--online-only)
+    #  ONLINE_ONLY=1
+    #  shift # past argument
+    #  #shift # past value
+    #  ;;
+    *)    # unknown option
+      #POSITIONAL+=("$1") # save it in an array for later
+      echo "unknown option $1"
+      print_the_help
+      shift # past argument
+      ;;
+  esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
+
+#if [[ -n "${LOG_ALL_PLOTS}" ]] ; then
+#  ONLINE_ONLY=
+#fi
+#echo "Number files in SEARCH PATH with EXTENSION:" $(ls -1 "${SEARCHPATH}"/*."${EXTENSION}" | wc -l)
+#if [[ -n $1 ]]; then
+#    echo "Last line of file specified as non-opt/last argument:"
+#    tail -1 "$1"
+#fi
+
+files=
+#if [[ "${CODATYPE}" == "coin" ]] ; then
+  raw_dir=/net/cdaq/cdaql3data/coda/data/raw
+  rawcopied_dir=/net/cdaq/cdaql3data/coda/data/raw.copiedtotape
+
+  if [[ -n "${USE_RAW_DIR}" ]] ; then
+    files=(${USE_RAW_DIR}/${CODATYPE}_*.dat)
+  else
+    files=(${rawcopied_dir}/${CODATYPE}_*.dat ${raw_dir}/${CODATYPE}_*.dat)
+  fi 
+
+  for afile in "${files[@]}"
+  do 
+    basename  $afile  .dat
+  done | sort  | tail -n 1 | sed "s/${CODATYPE}_all_//g" | sed 's/^0*//' 
+#fi
+
+

db2/auto_standard.kinematics

@@ -14146,3 +14146,63 @@ ppartmass = 0.1395706
 hpartmass = 0.0005109
 
 
+7105 - 7105
+gpbeam = 10.6
+gtargmass_amu = 1.00794
+htheta_lab = -21.09
+ptheta_lab = 13.11
+hpcentral = 0.0
+ppcentral = 5.925
+ppartmass = 0.1395706
+hpartmass = 0.0005109
+
+7105 - 7105
+gpbeam = 10.6
+gtargmass_amu = 1.00794
+htheta_lab = -21.09
+ptheta_lab = 13.11
+hpcentral = 0.0
+ppcentral = 5.925
+ppartmass = 0.1395706
+hpartmass = 0.0005109
+
+7105 - 7105
+gpbeam = 10.6
+gtargmass_amu = 1.00794
+htheta_lab = -21.09
+ptheta_lab = 13.11
+hpcentral = 0.0
+ppcentral = 5.925
+ppartmass = 0.1395706
+hpartmass = 0.0005109
+
+7105 - 7105
+gpbeam = 10.6
+gtargmass_amu = 1.00794
+htheta_lab = -21.09
+ptheta_lab = 13.11
+hpcentral = 0.0
+ppcentral = 5.925
+ppartmass = 0.1395706
+hpartmass = 0.0005109
+
+7105 - 7105
+gpbeam = 10.6
+gtargmass_amu = 1.00794
+htheta_lab = -21.09
+ptheta_lab = 13.1
+hpcentral = 0.0
+ppcentral = 5.925
+ppartmass = 0.1395706
+hpartmass = 0.0005109
+
+7105 - 7105
+gpbeam = 10.6
+gtargmass_amu = 1.00794
+htheta_lab = -21.09
+ptheta_lab = 13.1
+hpcentral = 0.0
+ppcentral = 5.925
+ppartmass = 0.1395706
+hpartmass = 0.0005109
+

db2/auto_standard_hms.kinematics

+2374 - 2374
+gpbeam = None
+gtargmass_amu = 12.0107
+htheta_lab = -21.0
+ptheta_lab = 25.0
+hpcentral = 4.1
+ppcentral = 3.5
+ppartmass = 0.0005109
+hpartmass = 0.0005109
+
+2376 - 2376
+gpbeam = None
+gtargmass_amu = 0
+htheta_lab = -21.0
+ptheta_lab = 25.05
+hpcentral = 4.1
+ppcentral = 3.5
+ppartmass = 0.0005109
+hpartmass = 0.0005109
+
+2376 - 2376
+gpbeam = 10.6
+gtargmass_amu = 1.00794
+htheta_lab = -21.09
+ptheta_lab = 13.11
+hpcentral = 0.0
+ppcentral = 5.925
+ppartmass = 0.0005109
+hpartmass = 0.0005109
+
+2376 - 2376
+gpbeam = 10.6
+gtargmass_amu = 12.0107
+htheta_lab = -21.09
+ptheta_lab = 13.1
+hpcentral = 0.0
+ppcentral = 5.925
+ppartmass = 0.0005109
+hpartmass = 0.0005109
+
+2376 - 2376
+gpbeam = 10.6
+gtargmass_amu = 12.0107
+htheta_lab = -21.09
+ptheta_lab = 13.1
+hpcentral = 0.0
+ppcentral = 5.925
+ppartmass = 0.0005109
+hpartmass = 0.0005109
+
+2376 - 2376
+gpbeam = 10.6
+gtargmass_amu = 12.0107
+htheta_lab = -21.09
+hpcentral = 0.0
+hpartmass = 0.0005109
+

online_bin/start_auto_replay

+#!/usr/bin/env python3
+
+from hallc.run_daemon import RunDaemon
+import subprocess
+import sys
+import os
+import libtmux
+
+p = None
+def launch_full_shms_replay(run_number):
+    global p
+    #command=' pwd && sleep 200 '
+    replay_dir = "/home/cdaq/whit/"
+    try:
+        replay_dir = os.environ['CURRENT_REPLAY_DIR']
+    except KeyError:
+        print("CURRENT_REPLAY_DIR is NOT SET! Why not?")
+        sys.exit(1)
+         
+    server = libtmux.Server(socket_name="replay")
+    sesh = server.find_where({"session_name":"full_replay"})
+    if not sesh:
+        sesh = server.new_session(session_name="full_replay",start_directory=replay_dir)
+ 
+    window = sesh.find_where({"window_name":"shms_run_{}".format(run_number)})
+    if window :
+        window.kill_window()
+    window = sesh.new_window(window_name="shms_run_{}".format(run_number), start_directory=replay_dir, attach=False,
+            window_shell="./bin/hc_shms_replay -r {} -n -1".format(run_number))
+    #p = subprocess.Popen(
+    #        [full_command ],
+    #        cwd='/home/cdaq/whit/online_jpsi/replay_jpsi/',
+    #        stdout = sys.stdout,
+    #        stderr = sys.stderr,
+    #        shell = True)
+        
+def launch_full_coin_replay(run_number):
+    global p
+    #command=' pwd && sleep 200 '
+    replay_dir = "/home/cdaq/whit/"
+    try:
+        replay_dir = os.environ['CURRENT_REPLAY_DIR']
+    except KeyError:
+        print("CURRENT_REPLAY_DIR is NOT SET! Why not?")
+        sys.exit(1)
+         
+    server = libtmux.Server(socket_name="replay")
+    sesh = server.find_where({"session_name":"full_replay"})
+    if not sesh:
+        sesh = server.new_session(session_name="full_replay",start_directory=replay_dir)
+ 
+    window = sesh.find_where({"window_name":"coin_run_{}".format(run_number)})
+    if window :
+        window.kill_window()
+    window = sesh.new_window(window_name="coin_run_{}".format(run_number), start_directory=replay_dir, attach=False,
+            window_shell="./bin/hc_coin_replay -r {} -n -1".format(run_number))
+
+def launch_full_hms_replay(run_number):
+    global p
+    #command=' pwd && sleep 200 '
+    replay_dir = "/home/cdaq/whit/"
+    try:
+        replay_dir = os.environ['CURRENT_REPLAY_DIR']
+    except KeyError:
+        print("CURRENT_REPLAY_DIR is NOT SET! Why not?")
+        sys.exit(1)
+         
+    server = libtmux.Server(socket_name="replay")
+    sesh = server.find_where({"session_name":"full_replay"})
+    if not sesh:
+        sesh = server.new_session(session_name="full_replay",start_directory=replay_dir)
+ 
+    window = sesh.find_where({"window_name":"hms_run_{}".format(run_number)})
+    if window :
+        window.kill_window()
+    window = sesh.new_window(window_name="hms_run_{}".format(run_number), start_directory=replay_dir, attach=False,
+            window_shell="./bin/hc_hms_replay -r {} -n -1".format(run_number))
+
+d = RunDaemon()
+d.on_event('run_start', launch_full_shms_replay, run_type='shms')
+d.on_event('run_start', launch_full_hms_replay, run_type='hms')
+d.on_event('run_start', launch_full_coin_replay, run_type='coin')
+d.start()
+

scripts/good_shms_counter.cxx

+#include "nlohmann/json.hpp"
+#include <cmath>
+#include <iostream>
+
+#include "ROOT/RDataFrame.hxx"
+#include "ROOT/RVec.hxx"
+
+#include "Math/Vector3D.h"
+#include "Math/Vector4D.h"
+#include "Math/VectorUtil.h"
+#include "TCanvas.h"
+#include "TLatex.h"
+#include "TStyle.h"
+#include "TSystem.h"
+R__LOAD_LIBRARY(libMathMore.so)
+R__LOAD_LIBRARY(libGenVector.so)
+
+#include "THStack.h"
+
+#ifdef __cpp_lib_filesystem
+#include <filesystem>
+namespace fs = std::filesystem;
+#else
+#include <experimental/filesystem>
+namespace fs = std::experimental::filesystem;
+#endif
+
+using Pvec3D = ROOT::Math::XYZVector;
+using Pvec4D = ROOT::Math::PxPyPzMVector;
+
+// VecOps::RVec is like std::vector with some extra bells and whistles
+using inters = ROOT::VecOps::RVec<int>;
+using doublers = ROOT::VecOps::RVec<double>;
+using floaters = ROOT::VecOps::RVec<float>;
+using shorters = ROOT::VecOps::RVec<short>;
+using nlohmann::json;
+
+void good_shms_counter(int RunNumber = 7146, int nevents = -1, int prompt =0, int update = 1) {
+
+  using nlohmann::json;
+  json j;
+  {
+    std::ifstream json_input_file("db2/run_list_shms.json");
+    try {
+      json_input_file >> j;
+    } catch(json::parse_error)  {
+      std::cerr << "error: json file, db2/run_list.json, is incomplete or has broken syntax.\n";
+      std::quick_exit(-127);
+    }
+  }
+
+  auto runnum_str = std::to_string(RunNumber);
+  if (j.find(runnum_str) == j.end()) {
+    std::cout << "Run " << RunNumber << " not found in db2/run_list_shms.json\n";
+    std::cout << "Check that run number and replay exists. \n";
+    std::cout << "If problem persists please contact Whit (717-341-1080)\n";
+  }
+  double P0_shms_setting =
+      j[runnum_str]["spectrometers"]["shms_momentum"].get<double>();
+  double P0_shms = std::abs(P0_shms_setting);
+
+  int ps2 = 0;
+  if (j[runnum_str].find("daq") != j[runnum_str].end()) {
+    ps2 = j[runnum_str]["daq"]["ps2"].get<int>();
+  } else {
+    std::cout << " using default ps2 = 0 \n";
+  }
+  //  The way the input rates are prescaled follows:
+  //       input-rate/(2^{val - 1} + 1)
+  double singles_ps_value = std::pow(2.0, ps2);
+  std::cout << "prescale value " << singles_ps_value << "\n";
+
+  std::string coda_type = "SHMS";
+  std::string rootfile = "ROOTfiles/";
+  rootfile += std::string("shms_replay_production_");
+  rootfile +=
+      std::to_string(RunNumber) + "_" + std::to_string(nevents) + ".root";
+
+  bool found_good_file = false; 
+  if(!gSystem->AccessPathName(rootfile.c_str())) {
+    TFile file(rootfile.c_str());
+    if (file.IsZombie()) {
+      std::cout << rootfile << " is a zombie.\n";
+      std::cout << " Did your replay finish?  Check that the it is done before running this script.\n";
+      //return;
+    } else {
+      found_good_file = true;
+    }
+  }
+  if(!found_good_file) {
+    rootfile = "ROOTfiles_online/";
+    rootfile += std::string("shms_replay_production_");
+    rootfile += std::to_string(RunNumber) + "_" + std::to_string(nevents) + ".root";
+
+    if(!gSystem->AccessPathName(rootfile.c_str())) {
+      TFile file(rootfile.c_str());
+      if (file.IsZombie()) {
+        std::cout << rootfile << " is a zombie.\n";
+        std::cout << " Did your replay finish?  Check that the it is done before running this script.\n";
+      } else {
+        found_good_file = true;
+      }
+    }
+  }
+  if(!found_good_file) {
+    std::cout << " Error: suitable root file not found\n";
+    return;
+  }
+
+  ROOT::EnableImplicitMT(24);
+
+  Pvec4D Pbeam(0, 0, 10.598, 0.000511);
+
+  //---------------------------------------------------------------------------
+  // Detector tree 
+  ROOT::RDataFrame d("T", rootfile);
+
+  // SHMS Scaler tree
+  ROOT::RDataFrame d_sh("TSP", rootfile);
+  //int N_scaler_events = *(d_sh.Count());
+
+  std::string hpdelta = "P.gtr.dp > -10 && P.gtr.dp < 10";
+  std::string epiCut = "P.ngcer.npeSum > 1.0 && P.cal.etottracknorm > 0.6 && "
+                       "P.cal.etottracknorm < 2.0 ";
+  // && H.cal.eprtracknorm  > 0.2
+  auto c_n_events_total = d.Count();
+  auto c_n_events_shms = d.Filter("(fEvtHdr.fEvtType == 6) || (fEvtHdr.fEvtType == 1)").Count();
+
+  auto d0 = d;
+  auto d0_coin_only = d.Filter("fEvtHdr.fEvtType == 6");
+  auto d0_shms_only  = d.Filter("fEvtHdr.fEvtType == 1");
+
+  // Apply the electron cuts
+  auto d_spec_cuts = d0.Filter(hpdelta).Filter(epiCut);
+  auto d_spec_cuts_coin_only = d_spec_cuts.Filter("fEvtHdr.fEvtType == 6");
+  auto d_spec_cuts_shms_only  = d_spec_cuts.Filter("fEvtHdr.fEvtType == 1");
+
+  //auto h_EOverP_0 = d0.Histo1D<doublers>(
+  //    {"shms_e_EoverP_0", "HMS total shower; SHMS E/P", 100, 0.05, 1.8}, "P.cal.etottracknorm");
+  //auto h_EOverP_2 = d_spec_cuts.Histo1D<doublers>(
+  //    {"hms_e_EoverP_2", "HMS total shower; SHMS E/P", 100, 0.05, 1.8}, "P.cal.etottracknorm");
+
+  //auto h_EprOverP_0 = d0.Histo1D<doublers>(
+  //    {"hms_e_EoverP_0", "HMS pre-shower; SHMS E/P", 100, 0.05, 1.8}, "P.cal.eprtracknorm");
+
+  //auto h_EprOverP_2 = d_spec_cuts.Histo1D<doublers>(
+  //    {"hms_e_EoverP_2", "HMS pre-shower; SHMS E/P", 100, 0.05, 1.8}, "P.cal.eprtracknorm");
+
+  auto h_event_type   = d0.Histo1D({"EvtType", "EvtType", 10, 0, 10}, "fEvtHdr.fEvtType");
+  auto h_event_type_2 = d_spec_cuts.Histo1D({"EvtType2", "EvtType", 10, 0, 10}, "fEvtHdr.fEvtType");
+
+  auto h_beta_0 = d_spec_cuts.Histo1D({"h_beta_0", "h_beta_0", 150, 0, 1.5}, "P.gtr.beta");
+
+  auto h_ngc_0 = d0         .Histo1D<doublers>( {"shms_ngcer", "SHMS NGC; nped", 100, 0, 30}, "P.ngcer.npeSum");
+  auto h_ngc_1 = d_spec_cuts.Histo1D<doublers>( {"shms_ngcer2", "SHMS NGC; nped", 100, 0, 30}, "P.ngcer.npeSum");
+
+  auto h_hgc_0 = d0         .Histo1D<doublers>( {"shms_hgcer", "SHMS HGC; nped", 100, 0, 30}, "P.hgcer.npeSum");
+  auto h_hgc_1 = d_spec_cuts.Histo1D<doublers>( {"shms_hgcer2","SHMS HGC; nped", 100, 0, 30}, "P.hgcer.npeSum");
+
+  auto bcm4b_charge        = d_sh.Max("P.BCM4B.scalerChargeCut");
+  auto el_real_scaler      = d_sh.Max("P.hEL_REAL.scaler");
+  auto time_1MHz           = d_sh.Max("P.1MHz.scalerTime");
+  auto time_1MHz_cut       = d_sh.Max("P.1MHz.scalerTimeCut");
+  //auto hTRIG1_ROC1_npassed = d_sh.Max("P.hTRIG1_ROC1.npassed");
+  //auto H_hTRIG1_scaler     = d_sh.Max("P.hTRIG1.scaler");
+
+  //{(hTRIG1_ROC1.npassed / H.hTRIG1.scaler)*100.0:%3.4f}
+  //H.hEL_REAL.scaler/P.1MHz.scalerTime)/1000
+  auto total_charge = bcm4b_charge;
+
+  auto c_e_yield_raw  = d.Count();
+  auto c_e_yield      = d_spec_cuts.Count();
+  auto c_T2_yield_raw = d0_shms_only.Count();
+  auto c_T6_yield_raw = d0_coin_only.Count();
+  auto c_T2_yield     = d_spec_cuts_shms_only.Count();
+  auto c_T6_yield     = d_spec_cuts_coin_only.Count();
+
+  // -------------------------------------
+  // End lazy eval
+  // -------------------------------------
+
+  double good_total_charge = *bcm4b_charge / 1000.0; // mC
+  //double shms_live_time = double(*hTRIG1_ROC1_npassed) / double(*H_hTRIG1_scaler);
+  double good_time = *time_1MHz_cut; // s
+
+  double shms_scaler_yield     = ((*el_real_scaler) / good_total_charge);
+  double shms_scaler_yield_unc = (std::sqrt(*el_real_scaler) / good_total_charge);
+
+  double shms_e_yield      = (*c_e_yield_raw) / (*total_charge);
+  double shms_e_yield2     = (*c_e_yield) / (*total_charge);
+
+  double ps_cor_shms_e_yield  = shms_e_yield *singles_ps_value;
+  double ps_cor_shms_e_yield2 = shms_e_yield2*singles_ps_value;
+
+  std::cout << " good_total_charge " << good_total_charge << " \n";
+  std::cout << " shms_scaler_yield " << shms_scaler_yield << " \n";
+
+  // Update counts list
+  json jruns;
+  {
+    std::ifstream input_file("db2/jpsi_run_count_list.json");
+    try {
+      input_file >> jruns;
+    } catch(json::parse_error)  {
+      std::cerr << "error: json file is incomplete or has broken syntax.\n";
+      std::quick_exit(-127);
+    }
+  }
+  std::string run_str                         = std::to_string(RunNumber);
+  jruns[run_str]["shms e raw counts"]         = int(*c_T2_yield_raw + *c_T6_yield_raw);
+  jruns[run_str]["shms e counts"]             = int(*c_T2_yield + *c_T6_yield);
+  jruns[run_str]["ps cor. shms e raw counts"] = int((*c_T2_yield_raw) * singles_ps_value + (*c_T6_yield_raw));
+  jruns[run_str]["ps cor. shms e counts"]     = int((*c_T2_yield_raw) * singles_ps_value + (*c_T6_yield_raw));
+  jruns[run_str]["charge bcm4b 2u cut"]       = good_total_charge;
+  jruns[run_str]["time 1MHz 2u cut"]          = good_time;
+  jruns[run_str]["shms e raw yield"]          = double((*c_T2_yield_raw) * singles_ps_value ) / good_total_charge;
+  jruns[run_str]["shms e yield"]              = double((*c_T2_yield) * singles_ps_value + (*c_T6_yield)) / good_total_charge;
+  jruns[run_str]["shms e yield unc."]         = double(std::sqrt(*c_T2_yield) * singles_ps_value + std::sqrt(*c_T6_yield)) / good_total_charge;
+  jruns[run_str]["shms ps2 factor"]           = singles_ps_value;
+  jruns[run_str]["shms scaler yield"]         = shms_scaler_yield;
+  jruns[run_str]["shms scaler yield unc."]    = shms_scaler_yield_unc;
+  //jruns[run_str]["hms live time"]         = hms_live_time;
+
+  if( update ) {
+    std::cout << "Updating db2/jpsi_run_count_list.json with shms counts\n" ;
+    std::ofstream json_output_file("db2/jpsi_run_count_list.json");
+    json_output_file << std::setw(4) << jruns << "\n";
+  }
+
+  std::cout << " ----------------------------------------------    \n";
+  std::cout << " SHMS Yield = "
+            << ps_cor_shms_e_yield2 << " = " << (*c_e_yield) << "/" << (*total_charge)
+            << "   counts/mC \n";
+  std::cout << " ----------------------------------------------    \n";
+  std::cout << " of  " << *c_n_events_total << " total triggers\n";
+  std::cout << " and " << *c_n_events_shms << " shms triggers\n";
+  // std::cout << " pions+kaons   : " << *coin_counts << "\n";
+  // std::cout << " pions         : " << *pion_count << "\n";
+  // std::cout << " random        : " << random_bg << "\n";
+  // std::cout << " ratio         : "  << pi_K_ratio << " \n";
+
+  // std::cout << " counts  : " << *hms_electron_counts << "\n";
+  // std::cout << " charge  : " << *total_charge << " uC\n";
+  // std::cout << " yield   : " << (*hms_electron_counts)  << " cnts, " <<
+  // shms_e_yield  << " cnts/uC\n"; std::cout << " singles : " <<
+  // (*hms_electron_counts2) << " cnts, " << shms_e_yield2 << " cnts/uC\n";
+  // std::cout << " coin    : " << (*coin_counts) << " cnts, " << coin_yield <<
+  // " cnts/uC\n";
+
+  // auto s_dc_x_fp           = d2.Histo1D({"s_dc_x_fp ","xy fp; x",100,-50,50},
+  // "P.dc.x_fp"); auto s_dc_y_fp           = d2.Histo1D({"s_dc_y_fp ","xy fp;
+  // y",100,-50,50}, "P.dc.y_fp"); auto s_dc_xp_fp          =
+  // d2.Histo1D({"s_dc_xp_fp","xy fp; xp",100,-50,50},"P.dc.xp_fp"); auto
+  // s_dc_yp_fp          = d2.Histo1D({"s_dc_yp_fp","xy fp;
+  // xp",100,-50,50},"P.dc.yp_fp");
+
+  // -----------------------------------------------------------
+  //
+  TCanvas *c = nullptr;
+  int b1 = 0;
+  int b2 = 0;
+  double hmax = 0.0;
+  THStack *hs = nullptr;
+  // TLatex latex;
+
+  gSystem->mkdir("results/good_shms_counter",true);
+
+  // ---------------------------------------------------------
+  //
+  c = new TCanvas();
+  c->Divide(2, 2);
+  c->cd(1);
+
+  // This call starts the loop over the data.
+  // A DrawCopy is used so that the histogram is not deleted at the end of
+  // scope, and thus stays visible on the canvas.
+  //h_EOverP_0->DrawCopy();
+
+  c->cd(2);
+  h_beta_0->DrawCopy();
+
+  c->cd(3);
+  gPad->SetLogy(true);
+  h_event_type->DrawCopy();
+  h_event_type_2->SetLineColor(2);
+  h_event_type_2->DrawCopy("same");
+
+  c->cd(4);
+  //h_coin_time->DrawCopy();
+  //h_coin_time2->SetLineColor(4);
+  //h_coin_time2->DrawCopy("same");
+  //h_coin_time3->SetLineColor(2);
+  //h_coin_time3->DrawCopy("same");
+
+  gPad->BuildLegend();
+
+  gSystem->mkdir("results/df_example", true);
+  c->SaveAs((std::string("results/df_example/c1_") + std::to_string(RunNumber) +
+             ".pdf")
+                .c_str());
+  c->SaveAs((std::string("results/df_example/c1_") + std::to_string(RunNumber) +
+             ".png")
+                .c_str());
+
+  // ---------------------------------------------------------
+  //
+  c = new TCanvas();
+  c->Divide(2, 2);
+  c->cd(1);
+  hs = new THStack("ngcer", "ngcer; npe");
+  hs->Add((TH1 *)h_ngc_0->Clone());
+  hs->Add((TH1 *)h_ngc_1->Clone());
+
+  hs->Draw("nostack");
+  gPad->SetLogy(true);
+  gPad->BuildLegend();
+
+  c->cd(2);
+  hs = new THStack("hgcer", "hgcer; npe");
+  hs->Add((TH1 *)h_hgc_0->Clone());
+  hs->Add((TH1 *)h_hgc_1->Clone());
+
+  hs->Draw("nostack");
+  gPad->SetLogy(true);
+  gPad->BuildLegend();
+
+
+
+  //// ---------------------------------------------------------
+  ////
+  //c = new TCanvas();
+  //hs = new THStack("SHMS_cal", "SHMS calorimeter; E/P");
+
+  //h_EOverP_0->SetLineColor(1);
+  //h_EOverP_2->SetLineColor(4);
+  //h_EprOverP_0->SetLineColor(8);
+  //h_EprOverP_2->SetLineColor(6);
+
+  //hs->Add((TH1 *)h_EOverP_0->Clone());
+  //hs->Add((TH1 *)h_EOverP_2->Clone());
+
+  //hs->Add((TH1 *)h_EprOverP_0->Clone());
+  //hs->Add((TH1 *)h_EprOverP_2->Clone());
+
+  //hs->Draw("nostack");
+  //gPad->SetLogy(true);
+  //gPad->BuildLegend();
+
+  //c->SaveAs((std::string("results/df_example/c2_") + std::to_string(RunNumber) +
+  //           ".pdf")
+  //              .c_str());
+  //c->SaveAs((std::string("results/df_example/c2_") + std::to_string(RunNumber) +
+  //           ".png")
+  //              .c_str());
+}

scripts/make_human_table.cxx

@@ -19,12 +19,12 @@ void make_human_table() {
   using nlohmann::json;
   json j;
   {
-    std::ifstream json_input_file("db2/run_list.json");
+    std::ifstream json_input_file("db2/run_list_shms.json");
     json_input_file >> j;
   }
   json j2;
   {
-    std::ifstream json_input_file("db2/run_count_list.json");
+    std::ifstream json_input_file("db2/jpsi_run_count_list.json");
     json_input_file >> j2;
   }
   json j3;

scripts/make_jpsi_table.cxx

+#include <fmt/core.h>
+#include <fmt/ostream.h>
+R__LOAD_LIBRARY(libfmt.so)
+
+#ifdef __cpp_lib_filesystem
+#include <filesystem>
+namespace fs = std::filesystem;
+#else
+#include <experimental/filesystem>
+namespace fs = std::experimental::filesystem;
+#endif
+
+#include "nlohmann/json.hpp"
+
+#include "TObject.h"
+
+void make_jpsi_table() {
+
+  using nlohmann::json;
+  json j;
+  {
+    std::ifstream json_input_file("db2/run_list_shms.json");
+    json_input_file >> j;
+  }
+  json j2;
+  {
+    std::ifstream json_input_file("db2/jpsi_run_count_list.json");
+    json_input_file >> j2;
+  }
+  json j3;
+  {
+    std::ifstream json_input_file("db2/run_list_extra.json");
+    try {
+      json_input_file >> j3;
+    } catch (json::parse_error) {
+      std::cerr << "error: json file, db2/run_list.json, is incomplete or has broken syntax.\n";
+      std::quick_exit(-127);
+    }
+  }
+  std::ofstream runs_text("jpsi_runs.txt");
+
+  // std::cout << j.dump(2);
+  auto print_header = []() {
+    // print header
+    std::cout << "\n";
+    fmt::print(" {:<5} ", "Run");
+    fmt::print(" {:^5} ", "Target");
+    fmt::print(" {:^7} ", "P_hms");
+    fmt::print(" {:^7} ", "th_hms");
+    fmt::print(" {:^7} ", "P_shms");
+    fmt::print(" {:^7} ", "th_shms");
+    fmt::print(" {:^21} ", "start time");
+    fmt::print(" {:^21} ", "end time");
+    fmt::print(" {:>9} ", "count");
+    fmt::print(" {:>9} ", "charge");
+    fmt::print(" {:>9} ", "yield");
+    fmt::print(" {:>9} ", "triggers");
+    fmt::print(" {:<} ", "comment");
+    std::cout << "\n";
+  };
+
+  std::cout << " runs : \n";
+
+  std::string old_target = "";
+
+  double p_hms   = 0.0;
+  double th_hms  = 0.0;
+  double p_shms  = 0.0;
+  double th_shms = 0.0;
+
+  for (json::iterator it = j.begin(); it != j.end(); ++it) {
+    auto runjs = it.value();
+
+    std::string target_lab = runjs["target"]["target_label"].get<std::string>();
+    if (target_lab != old_target) {
+      print_header();
+    }
+
+    if ((p_hms != runjs["spectrometers"]["hms_momentum"].get<double>()) ||
+        (th_hms != runjs["spectrometers"]["hms_angle"].get<double>()) ||
+        (p_shms != runjs["spectrometers"]["shms_momentum"].get<double>()) ||
+        (th_shms != runjs["spectrometers"]["shms_angle"].get<double>())) {
+      print_header();
+    }
+
+    p_hms   = runjs["spectrometers"]["hms_momentum"].get<double>();
+    th_hms  = runjs["spectrometers"]["hms_angle"].get<double>();
+    p_shms  = runjs["spectrometers"]["shms_momentum"].get<double>();
+    th_shms = runjs["spectrometers"]["shms_angle"].get<double>();
+
+    old_target = target_lab;
+
+    runs_text << std::stoi(it.key()) << "\n";
+
+    fmt::print(" {:<5} ", std::stoi(it.key()));
+    fmt::print(" {:^5} ", target_lab);
+    fmt::print(" {:>7.3f} ", runjs["spectrometers"]["hms_momentum"].get<double>());
+    fmt::print(" {:>7.2f} ", runjs["spectrometers"]["hms_angle"].get<double>());
+    fmt::print(" {:>7.3f} ", runjs["spectrometers"]["shms_momentum"].get<double>());
+    fmt::print(" {:>7.2f} ", runjs["spectrometers"]["shms_angle"].get<double>());
+    if (runjs["run_info"].find("start_time") != runjs["run_info"].end()) {
+      fmt::print(" {:^21} ", runjs["run_info"]["start_time"].get<std::string>());
+    } else {
+      fmt::print(" {:21} ", "");
+    }
+    if (runjs["run_info"].find("end_time") != runjs["run_info"].end()) {
+      fmt::print(" {:^21} ", runjs["run_info"]["end_time"].get<std::string>());
+    } else {
+      fmt::print(" {:21} ", "");
+    }
+
+    double total_charge = 0.0001;
+    if (runjs.find("total_charge") != runjs.end()) {
+      total_charge = runjs["total_charge"].get<double>() / 1000.0;
+      // fmt::print("/ {:>10.1f} ", total_charge);
+    }
+    // else {
+    //  fmt::print(" {:>11} ", "");
+    //}
+    if (j2.count(it.key()) != 0) {
+      try {
+
+        double shms_yield = j2[it.key()]["shms e counts"].get<double>();
+        // double pi_yield = j2[it.key()]["pion bg sub. counts"].get<double>();
+        fmt::print(" {:>9.1f} ", shms_yield);
+        if (j2[it.key()].find("charge bcm4b 2u cut") != j2[it.key()].end()) {
+          total_charge = j2[it.key()]["charge bcm4b 2u cut"].get<double>();
+          // total_charge = runjs["total_charge"].get<double>() / 1000.0;
+          fmt::print("/ {:<6.1f} ", total_charge);
+        } else {
+          fmt::print(" {:>11} ", "");
+        }
+        fmt::print(" {:>9.1f} ", shms_yield / total_charge);
+        int n_events = j2[it.key()]["total trigger events"].get<int>();
+        fmt::print(" {:>9d} ", n_events);
+      } catch (std::domain_error) {
+        // you suck
+        //std::cout << "you suck\n";
+      }
+      // double live_time = j2[it.key()]["live time"].get<int>();
+    } else {
+      fmt::print(" {:>9} ", "");
+      fmt::print(" {:>9} ", "");
+      fmt::print(" {:>10} ", "");
+    }
+    std::string comment;
+    if (j3.count(it.key()) != 0) {
+      if (j3[it.key()].find("comment") != j3[it.key()].end()) {
+        try {
+          comment = j3[it.key()]["comment"].get<std::string>();
+        } catch (std::domain_error) {
+          // you suck
+        }
+      }
+    }
+    fmt::print(" {:<} ", comment);
+    std::cout << "\n";
+  }
+  print_header();
+}

scripts/muon_pion_separation.cxx

+#include "nlohmann/json.hpp"
+#include <cmath>
+#include <iostream>
+
+#include "ROOT/RDataFrame.hxx"
+#include "ROOT/RVec.hxx"
+
+#include "Math/Vector3D.h"
+#include "Math/Vector4D.h"
+#include "Math/VectorUtil.h"
+#include "TCanvas.h"
+#include "TLatex.h"
+#include "TStyle.h"
+#include "TSystem.h"
+R__LOAD_LIBRARY(libMathMore.so)
+R__LOAD_LIBRARY(libGenVector.so)
+
+#include "THStack.h"
+
+#ifdef __cpp_lib_filesystem
+#include <filesystem>
+namespace fs = std::filesystem;
+#else
+#include <experimental/filesystem>
+namespace fs = std::experimental::filesystem;
+#endif
+
+using Pvec3D = ROOT::Math::XYZVector;
+using Pvec4D = ROOT::Math::PxPyPzMVector;
+
+// VecOps::RVec is like std::vector with some extra bells and whistles
+using inters = ROOT::VecOps::RVec<int>;
+using doublers = ROOT::VecOps::RVec<double>;
+using floaters = ROOT::VecOps::RVec<float>;
+using shorters = ROOT::VecOps::RVec<short>;
+using nlohmann::json;
+
+void muon_pion_separation(int RunNumber = 7146, int nevents = -1, int prompt =0, int update = 1) {
+
+  using nlohmann::json;
+  json j;
+  {
+    std::ifstream json_input_file("db2/run_list_shms.json");
+    try {
+      json_input_file >> j;
+    } catch(json::parse_error)  {
+      std::cerr << "error: json file, db2/run_list.json, is incomplete or has broken syntax.\n";
+      std::quick_exit(-127);
+    }
+  }
+
+  auto runnum_str = std::to_string(RunNumber);
+  if (j.find(runnum_str) == j.end()) {
+    std::cout << "Run " << RunNumber << " not found in db2/run_list_shms.json\n";
+    std::cout << "Check that run number and replay exists. \n";
+    std::cout << "If problem persists please contact Whit (717-341-1080)\n";
+  }
+  double P0_shms_setting =
+      j[runnum_str]["spectrometers"]["shms_momentum"].get<double>();
+  double P0_shms = std::abs(P0_shms_setting);
+
+  int ps2 = 0;
+  if (j[runnum_str].find("daq") != j[runnum_str].end()) {
+    ps2 = j[runnum_str]["daq"]["ps2"].get<int>();
+  } else {
+    std::cout << " using default ps2 = 0 \n";
+  }
+  //  The way the input rates are prescaled follows:
+  //       input-rate/(2^{val - 1} + 1)
+  double singles_ps_value = std::pow(2.0, ps2);
+  std::cout << "prescale value " << singles_ps_value << "\n";
+
+  std::string coda_type = "SHMS";
+  std::string rootfile = "ROOTfiles/";
+  rootfile += std::string("shms_replay_production_");
+  rootfile +=
+      std::to_string(RunNumber) + "_" + std::to_string(nevents) + ".root";
+
+  bool found_good_file = false; 
+  if(!gSystem->AccessPathName(rootfile.c_str())) {
+    TFile file(rootfile.c_str());
+    if (file.IsZombie()) {
+      std::cout << rootfile << " is a zombie.\n";
+      std::cout << " Did your replay finish?  Check that the it is done before running this script.\n";
+      //return;
+    } else {
+      found_good_file = true;
+    }
+  }
+  if(!found_good_file) {
+    rootfile = "ROOTfiles_online/";
+    rootfile += std::string("shms_replay_production_");
+    rootfile += std::to_string(RunNumber) + "_" + std::to_string(nevents) + ".root";
+
+    if(!gSystem->AccessPathName(rootfile.c_str())) {
+      TFile file(rootfile.c_str());
+      if (file.IsZombie()) {
+        std::cout << rootfile << " is a zombie.\n";
+        std::cout << " Did your replay finish?  Check that the it is done before running this script.\n";
+      } else {
+        found_good_file = true;
+      }
+    }
+  }
+  if(!found_good_file) {
+    std::cout << " Error: suitable root file not found\n";
+    return;
+  }
+
+  ROOT::EnableImplicitMT(24);
+
+  Pvec4D Pbeam(0, 0, 10.598, 0.000511);
+
+  //---------------------------------------------------------------------------
+  // Detector tree 
+  ROOT::RDataFrame d("T", rootfile);
+
+  // SHMS Scaler tree
+  ROOT::RDataFrame d_sh("TSP", rootfile);
+  //int N_scaler_events = *(d_sh.Count());
+
+  std::string hpdelta = "P.gtr.dp > -10 && P.gtr.dp < 10";
+  std::string epiCut = "P.ngcer.npeSum > 1.0 && P.cal.etottracknorm > 0.6 && "
+                       "P.cal.etottracknorm < 2.0 ";
+  // && H.cal.eprtracknorm  > 0.2
+  auto c_n_events_total = d.Count();
+  auto c_n_events_shms = d.Filter("(fEvtHdr.fEvtType == 6) || (fEvtHdr.fEvtType == 1)").Count();
+
+  auto d0 = d;
+  auto d0_coin_only = d.Filter("fEvtHdr.fEvtType == 6");
+  auto d0_shms_only  = d.Filter("fEvtHdr.fEvtType == 1");
+
+  // Apply the electron cuts
+  auto d_spec_cuts = d0.Filter(hpdelta).Filter(epiCut);
+  auto d_spec_cuts_coin_only = d_spec_cuts.Filter("fEvtHdr.fEvtType == 6");
+  auto d_spec_cuts_shms_only  = d_spec_cuts.Filter("fEvtHdr.fEvtType == 1");
+
+  //auto h_EOverP_0 = d0.Histo1D<doublers>(
+  //    {"shms_e_EoverP_0", "HMS total shower; SHMS E/P", 100, 0.05, 1.8}, "P.cal.etottracknorm");
+  //auto h_EOverP_2 = d_spec_cuts.Histo1D<doublers>(
+  //    {"hms_e_EoverP_2", "HMS total shower; SHMS E/P", 100, 0.05, 1.8}, "P.cal.etottracknorm");
+
+  //auto h_EprOverP_0 = d0.Histo1D<doublers>(
+  //    {"hms_e_EoverP_0", "HMS pre-shower; SHMS E/P", 100, 0.05, 1.8}, "P.cal.eprtracknorm");
+
+  //auto h_EprOverP_2 = d_spec_cuts.Histo1D<doublers>(
+  //    {"hms_e_EoverP_2", "HMS pre-shower; SHMS E/P", 100, 0.05, 1.8}, "P.cal.eprtracknorm");
+
+  auto h_event_type   = d0.Histo1D({"EvtType", "EvtType", 10, 0, 10}, "fEvtHdr.fEvtType");
+  auto h_event_type_2 = d_spec_cuts.Histo1D({"EvtType2", "EvtType", 10, 0, 10}, "fEvtHdr.fEvtType");
+
+  auto h_beta_0 = d_spec_cuts.Histo1D({"h_beta_0", "h_beta_0", 150, 0, 1.5}, "P.gtr.beta");
+
+  auto h_ngc_0 = d0         .Histo1D<doublers>( {"shms_ngcer", "SHMS NGC; nped", 100, 0, 30}, "P.ngcer.npeSum");
+  auto h_ngc_1 = d_spec_cuts.Histo1D<doublers>( {"shms_ngcer2", "SHMS NGC; nped", 100, 0, 30}, "P.ngcer.npeSum");
+
+  auto h_hgc_0 = d0         .Histo1D<doublers>( {"shms_hgcer", "SHMS HGC; nped", 100, 0, 30}, "P.hgcer.npeSum");
+  auto h_hgc_1 = d_spec_cuts.Histo1D<doublers>( {"shms_hgcer2","SHMS HGC; nped", 100, 0, 30}, "P.hgcer.npeSum");
+
+  auto bcm4b_charge        = d_sh.Max("P.BCM4B.scalerChargeCut");
+  auto el_real_scaler      = d_sh.Max("P.hEL_REAL.scaler");
+  auto time_1MHz           = d_sh.Max("P.1MHz.scalerTime");
+  auto time_1MHz_cut       = d_sh.Max("P.1MHz.scalerTimeCut");
+  //auto hTRIG1_ROC1_npassed = d_sh.Max("P.hTRIG1_ROC1.npassed");
+  //auto H_hTRIG1_scaler     = d_sh.Max("P.hTRIG1.scaler");
+
+  //{(hTRIG1_ROC1.npassed / H.hTRIG1.scaler)*100.0:%3.4f}
+  //H.hEL_REAL.scaler/P.1MHz.scalerTime)/1000
+  auto total_charge = bcm4b_charge;
+
+  auto c_e_yield_raw  = d.Count();
+  auto c_e_yield      = d_spec_cuts.Count();
+  auto c_T2_yield_raw = d0_shms_only.Count();
+  auto c_T6_yield_raw = d0_coin_only.Count();
+  auto c_T2_yield     = d_spec_cuts_shms_only.Count();
+  auto c_T6_yield     = d_spec_cuts_coin_only.Count();
+
+  // -------------------------------------
+  // End lazy eval
+  // -------------------------------------
+
+  double good_total_charge = *bcm4b_charge / 1000.0; // mC
+  //double shms_live_time = double(*hTRIG1_ROC1_npassed) / double(*H_hTRIG1_scaler);
+  double good_time = *time_1MHz_cut; // s
+
+  double shms_scaler_yield     = ((*el_real_scaler) / good_total_charge);
+  double shms_scaler_yield_unc = (std::sqrt(*el_real_scaler) / good_total_charge);
+
+  double shms_e_yield      = (*c_e_yield_raw) / (*total_charge);
+  double shms_e_yield2     = (*c_e_yield) / (*total_charge);
+
+  double ps_cor_shms_e_yield  = shms_e_yield *singles_ps_value;
+  double ps_cor_shms_e_yield2 = shms_e_yield2*singles_ps_value;
+
+  std::cout << " good_total_charge " << good_total_charge << " \n";
+  std::cout << " shms_scaler_yield " << shms_scaler_yield << " \n";
+
+  // Update counts list
+  json jruns;
+  {
+    std::ifstream input_file("db2/jpsi_run_count_list.json");
+    try {
+      input_file >> jruns;
+    } catch(json::parse_error)  {
+      std::cerr << "error: json file is incomplete or has broken syntax.\n";
+      std::quick_exit(-127);
+    }
+  }
+  std::string run_str                         = std::to_string(RunNumber);
+  jruns[run_str]["shms e raw counts"]         = int(*c_T2_yield_raw + *c_T6_yield_raw);
+  jruns[run_str]["shms e counts"]             = int(*c_T2_yield + *c_T6_yield);
+  jruns[run_str]["ps cor. shms e raw counts"] = int((*c_T2_yield_raw) * singles_ps_value + (*c_T6_yield_raw));
+  jruns[run_str]["ps cor. shms e counts"]     = int((*c_T2_yield_raw) * singles_ps_value + (*c_T6_yield_raw));
+  jruns[run_str]["charge bcm4b 2u cut"]       = good_total_charge;
+  jruns[run_str]["time 1MHz 2u cut"]          = good_time;
+  jruns[run_str]["shms e raw yield"]          = double((*c_T2_yield_raw) * singles_ps_value ) / good_total_charge;
+  jruns[run_str]["shms e yield"]              = double((*c_T2_yield) * singles_ps_value + (*c_T6_yield)) / good_total_charge;
+  jruns[run_str]["shms e yield unc."]         = double(std::sqrt(*c_T2_yield) * singles_ps_value + std::sqrt(*c_T6_yield)) / good_total_charge;
+  jruns[run_str]["shms ps2 factor"]           = singles_ps_value;
+  jruns[run_str]["shms scaler yield"]         = shms_scaler_yield;
+  jruns[run_str]["shms scaler yield unc."]    = shms_scaler_yield_unc;
+  //jruns[run_str]["hms live time"]         = hms_live_time;
+
+  if( update ) {
+    std::cout << "Updating db2/jpsi_run_count_list.json with shms counts\n" ;
+    std::ofstream json_output_file("db2/jpsi_run_count_list.json");
+    json_output_file << std::setw(4) << jruns << "\n";
+  }
+
+  std::cout << " ----------------------------------------------    \n";
+  std::cout << " SHMS Yield = "
+            << ps_cor_shms_e_yield2 << " = " << (*c_e_yield) << "/" << (*total_charge)
+            << "   counts/mC \n";
+  std::cout << " ----------------------------------------------    \n";
+  std::cout << " of  " << *c_n_events_total << " total triggers\n";
+  std::cout << " and " << *c_n_events_shms << " shms triggers\n";
+  // std::cout << " pions+kaons   : " << *coin_counts << "\n";
+  // std::cout << " pions         : " << *pion_count << "\n";
+  // std::cout << " random        : " << random_bg << "\n";
+  // std::cout << " ratio         : "  << pi_K_ratio << " \n";
+
+  // std::cout << " counts  : " << *hms_electron_counts << "\n";
+  // std::cout << " charge  : " << *total_charge << " uC\n";
+  // std::cout << " yield   : " << (*hms_electron_counts)  << " cnts, " <<
+  // shms_e_yield  << " cnts/uC\n"; std::cout << " singles : " <<
+  // (*hms_electron_counts2) << " cnts, " << shms_e_yield2 << " cnts/uC\n";
+  // std::cout << " coin    : " << (*coin_counts) << " cnts, " << coin_yield <<
+  // " cnts/uC\n";
+
+  // auto s_dc_x_fp           = d2.Histo1D({"s_dc_x_fp ","xy fp; x",100,-50,50},
+  // "P.dc.x_fp"); auto s_dc_y_fp           = d2.Histo1D({"s_dc_y_fp ","xy fp;
+  // y",100,-50,50}, "P.dc.y_fp"); auto s_dc_xp_fp          =
+  // d2.Histo1D({"s_dc_xp_fp","xy fp; xp",100,-50,50},"P.dc.xp_fp"); auto
+  // s_dc_yp_fp          = d2.Histo1D({"s_dc_yp_fp","xy fp;
+  // xp",100,-50,50},"P.dc.yp_fp");
+
+  // -----------------------------------------------------------
+  //
+  TCanvas *c = nullptr;
+  int b1 = 0;
+  int b2 = 0;
+  double hmax = 0.0;
+  THStack *hs = nullptr;
+  // TLatex latex;
+
+  gSystem->mkdir("results/good_shms_counter",true);
+
+  // ---------------------------------------------------------
+  //
+  c = new TCanvas();
+  c->Divide(2, 2);
+  c->cd(1);
+
+  // This call starts the loop over the data.
+  // A DrawCopy is used so that the histogram is not deleted at the end of
+  // scope, and thus stays visible on the canvas.
+  //h_EOverP_0->DrawCopy();
+
+  c->cd(2);
+  h_beta_0->DrawCopy();
+
+  c->cd(3);
+  gPad->SetLogy(true);
+  h_event_type->DrawCopy();
+  h_event_type_2->SetLineColor(2);
+  h_event_type_2->DrawCopy("same");
+
+  c->cd(4);
+  //h_coin_time->DrawCopy();
+  //h_coin_time2->SetLineColor(4);
+  //h_coin_time2->DrawCopy("same");
+  //h_coin_time3->SetLineColor(2);
+  //h_coin_time3->DrawCopy("same");
+
+  gPad->BuildLegend();
+
+  gSystem->mkdir("results/df_example", true);
+  c->SaveAs((std::string("results/df_example/c1_") + std::to_string(RunNumber) +
+             ".pdf")
+                .c_str());
+  c->SaveAs((std::string("results/df_example/c1_") + std::to_string(RunNumber) +
+             ".png")
+                .c_str());
+
+  // ---------------------------------------------------------
+  //
+  c = new TCanvas();
+  c->Divide(2, 2);
+  c->cd(1);
+  hs = new THStack("ngcer", "ngcer; npe");
+  hs->Add((TH1 *)h_ngc_0->Clone());
+  hs->Add((TH1 *)h_ngc_1->Clone());
+
+  hs->Draw("nostack");
+  gPad->SetLogy(true);
+  gPad->BuildLegend();
+
+  c->cd(2);
+  hs = new THStack("hgcer", "hgcer; npe");
+  hs->Add((TH1 *)h_hgc_0->Clone());
+  hs->Add((TH1 *)h_hgc_1->Clone());
+
+  hs->Draw("nostack");
+  gPad->SetLogy(true);
+  gPad->BuildLegend();
+
+
+
+  //// ---------------------------------------------------------
+  ////
+  //c = new TCanvas();
+  //hs = new THStack("SHMS_cal", "SHMS calorimeter; E/P");
+
+  //h_EOverP_0->SetLineColor(1);
+  //h_EOverP_2->SetLineColor(4);
+  //h_EprOverP_0->SetLineColor(8);
+  //h_EprOverP_2->SetLineColor(6);
+
+  //hs->Add((TH1 *)h_EOverP_0->Clone());
+  //hs->Add((TH1 *)h_EOverP_2->Clone());
+
+  //hs->Add((TH1 *)h_EprOverP_0->Clone());
+  //hs->Add((TH1 *)h_EprOverP_2->Clone());
+
+  //hs->Draw("nostack");
+  //gPad->SetLogy(true);
+  //gPad->BuildLegend();
+
+  //c->SaveAs((std::string("results/df_example/c2_") + std::to_string(RunNumber) +
+  //           ".pdf")
+  //              .c_str());
+  //c->SaveAs((std::string("results/df_example/c2_") + std::to_string(RunNumber) +
+  //           ".png")
+  //              .c_str());
+}

scripts/plot_shms_singles.cxx

+#ifdef __cpp_lib_filesystem
+#include <filesystem>
+namespace fs = std::filesystem;
+#else
+#include <experimental/filesystem>
+namespace fs = std::experimental::filesystem;
+#endif
+
+#include "nlohmann/json.hpp"
+
+#include "TObject.h"
+#include "TGraph.h"
+
+
+void plot_shms_singles(int start_run = 0){
+
+  using nlohmann::json;
+  json j;
+  {
+    std::ifstream json_input_file("db2/jpsi_run_count_list.json");
+    try {
+      json_input_file >> j;
+    } catch(json::parse_error) {
+      std::cerr << "error: json file, db2/run_count_list.json, is incomplete or has broken syntax.\n";
+      std::quick_exit(-127);
+    }
+  }
+  json j2;
+  {
+    std::ifstream json_input_file("db2/run_list_shms.json");
+    try {
+      json_input_file >> j2;
+    } catch(json::parse_error)  {
+      std::cerr << "error: json file, db2/run_list.json, is incomplete or has broken syntax.\n";
+      std::quick_exit(-127);
+    }
+  }
+
+
+  std::vector<double> runs;
+  std::vector<double> hms_yields;
+  std::vector<double> hms_yield_uncs;
+
+  std::vector<double> runs2;
+  std::vector<double> hms_yields2;
+  std::vector<double> hms_yield_uncs2;
+
+  std::vector<double> runs3;
+  std::vector<double> yield_ratios;
+  std::vector<double> yield_ratio_uncs;
+
+  std::vector<double> target_ids;
+
+  TGraphErrors* gr_hms_count  = new TGraphErrors();
+  TGraphErrors* gr_hms_scaler = new TGraphErrors();
+  TGraphErrors* gr_yield_ratio = new TGraphErrors();
+
+
+  double scaler_reduce_factor = 0.4;
+
+  for (json::iterator it = j.begin(); it != j.end(); ++it) {
+
+    if (j2.find(it.key()) == j2.end()) {
+      continue;
+    }
+    if (j2[it.key()]["target"]["target_id"].get<int>() != 3) {
+      continue;
+    }
+
+    int run_number = std::stoi(it.key());
+    int run_number2 = 0;
+    double yield      = 0.0;
+    double yield2      =0.0;
+
+    if(run_number < start_run ) {
+      continue;
+    }
+
+    if (it.value().find("shms e yield") != it.value().end()) {
+      try {
+      yield      = it.value()["shms e yield"].get<double>();
+      runs.push_back(double(run_number));
+      hms_yields.push_back(yield);
+      if (it.value().find("shms e yield unc.") != it.value().end()) {
+        hms_yield_uncs.push_back(it.value()["shms e yield unc."].get<double>());
+      } else {
+        hms_yield_uncs.push_back(0.0);
+      }
+      } catch(std::domain_error ) {
+        continue;
+        //you suck
+      }
+    }
+    if (it.value().find("shms scaler yield") != it.value().end()) {
+      run_number2 = std::stoi(it.key());
+      yield2      = it.value()["shms scaler yield"].get<double>();
+      runs2.push_back(double(run_number2));
+      hms_yields2.push_back(scaler_reduce_factor*yield2);
+      if (it.value().find("shms scaler yield unc.") != it.value().end()) {
+        hms_yield_uncs2.push_back(scaler_reduce_factor*it.value()["shms scaler yield unc."].get<double>());
+      } else {
+        hms_yield_uncs2.push_back(0.0);
+      }
+    }
+    if (it.value().find("shms scaler yield") != it.value().end()) {
+      if (it.value().find("shms e yield") != it.value().end()) {
+        runs3.push_back(double(run_number2));
+        double yield_ratio = yield / yield2;
+        yield_ratios.push_back(yield_ratio);
+        if ((it.value().find("shms e yield unc.") != it.value().end()) &&
+            (it.value().find("shms scaler yield unc.") != it.value().end())) {
+          double y1 = it.value()["shms e yield unc."].get<double>();
+          double y2 = it.value()["shms scaler yield unc."].get<double>();
+          yield_ratio_uncs.push_back(yield_ratio*std::sqrt(y1 * y1/(yield*yield) + y2 * y2/(yield2*yield2)));
+        } else {
+          yield_ratio_uncs.push_back(0.0);
+        }
+      }
+    }
+  }
+
+  // ------------------------------
+  //
+  TGraphErrors* gr = nullptr;
+  TMultiGraph*  mg = nullptr;
+  TCanvas*      c  = nullptr;
+
+  // ------------------------------
+  //
+  c  = new TCanvas("c1","c1",1600,1200);
+  c->Divide(1, 2);
+
+  c->cd(1);
+  mg = new TMultiGraph();
+  mg->SetTitle("; run number ; counts/mC");
+
+  std::vector<double> zeros(runs.size());
+  gr = new TGraphErrors(runs.size(), runs.data(), hms_yields.data(),zeros.data(),hms_yield_uncs.data());
+  gr->SetMarkerStyle(20);
+  mg->Add(gr,"p");
+
+  std::vector<double> zeros2(runs2.size());
+  gr = new TGraphErrors(runs2.size(), runs2.data(), hms_yields2.data(),zeros2.data(),hms_yield_uncs2.data());
+  gr->SetMarkerStyle(23);
+  gr->SetMarkerColor(2);
+  mg->Add(gr,"p");
+
+  mg->Draw("a");
+
+  c->cd(2);
+  std::vector<double> zeros3(runs3.size());
+  gr = new TGraphErrors(runs3.size(), runs3.data(), yield_ratios.data(),zeros3.data(),yield_ratio_uncs.data());
+  gr->SetTitle("; run number ; ps4*T2 + T6 counts / scaler");
+  gr->SetMarkerStyle(20);
+  //gr->Draw("alp");
+  //gr->GetYaxis()->SetRangeUser(0.3, 0.45);
+  gr->Draw("alp");
+
+  c->SaveAs("results/hms_singles_vs_run_number.png");
+  c->SaveAs("results/hms_singles_vs_run_number.pdf");
+}

scripts/replay_production_shms.cxx

+// ----------------------------------------------
+R__LOAD_LIBRARY(libsimple_epics.so)
+#include "simple_epics/PVList.h"
+
+R__LOAD_LIBRARY(libScandalizer.so)
+#include "scandalizer/PostProcessors.h"
+
+void replay_production_shms(Int_t RunNumber = 0, Int_t MaxEvent = 0) {
+
+  spdlog::set_level(spdlog::level::warn);
+  spdlog::flush_every(std::chrono::seconds(5));
+
+
+  // Get RunNumber and MaxEvent if not provided.
+  if(RunNumber == 0) {
+    cout << "Enter a Run Number (-1 to exit): ";
+    if( RunNumber<=0 ) return;
+  }
+  if(MaxEvent == 0) {
+    cout << "\nNumber of Events to analyze: ";
+    cin >> MaxEvent;
+    if(MaxEvent == 0) {
+      cerr << "...Invalid entry\n";
+      exit;
+    }
+  }
+
+  // Create file name patterns.
+  const char* RunFileNamePattern = "shms_all_%05d.dat";
+  vector<TString> pathList;
+  pathList.push_back(".");
+  pathList.push_back("./raw_coda");
+  pathList.push_back("./raw");
+  pathList.push_back("./raw.copiedtotape");
+  pathList.push_back("./raw/../raw.copiedtotape");
+  pathList.push_back("./cache");
+
+  //const char* RunFileNamePattern = "raw/coin_all_%05d.dat";
+  const char* ROOTFileNamePattern = "ROOTfiles/shms_replay_production_%d_%d.root";
+
+  // Load global parameters
+  gHcParms->Define("gen_run_number", "Run Number", RunNumber);
+  gHcParms->AddString("g_ctp_database_filename", "DBASE/SHMS/standard.database");
+  gHcParms->Load(gHcParms->GetString("g_ctp_database_filename"), RunNumber);
+  gHcParms->Load(gHcParms->GetString("g_ctp_parm_filename"));
+  gHcParms->Load(gHcParms->GetString("g_ctp_kinematics_filename"), RunNumber);
+  // Load parameters for SHMS trigger configuration
+  gHcParms->Load("PARAM/TRIG/tshms.param");
+  // Load fadc debug parameters
+  gHcParms->Load("PARAM/SHMS/GEN/p_fadc_debug.param");
+
+  // Load the Hall C detector map
+  gHcDetectorMap = new THcDetectorMap();
+  gHcDetectorMap->Load("MAPS/SHMS/DETEC/STACK/shms_stack.map");
+  
+  // ========
+  //  SHMS 
+  // ========
+
+  // Set up the equipment to be analyzed.
+  THcHallCSpectrometer* SHMS = new THcHallCSpectrometer("P", "SHMS");
+  SHMS->SetEvtType(1);
+  SHMS->AddEvtType(4);
+  SHMS->AddEvtType(5);
+  SHMS->AddEvtType(6);
+  SHMS->AddEvtType(7);
+  gHaApps->Add(SHMS);
+  // Add Noble Gas Cherenkov to SHMS apparatus
+  THcCherenkov* pngcer = new THcCherenkov("ngcer", "Noble Gas Cherenkov");
+  SHMS->AddDetector(pngcer);
+  // Add drift chambers to SHMS apparatus
+  THcDC* pdc = new THcDC("dc", "Drift Chambers");
+  SHMS->AddDetector(pdc);
+  // Add hodoscope to SHMS apparatus
+  THcHodoscope* phod = new THcHodoscope("hod", "Hodoscope");
+  SHMS->AddDetector(phod);
+  // Add Heavy Gas Cherenkov to SHMS apparatus
+  THcCherenkov* phgcer = new THcCherenkov("hgcer", "Heavy Gas Cherenkov");
+  SHMS->AddDetector(phgcer);
+  // Add Aerogel Cherenkov to SHMS apparatus
+  THcAerogel* paero = new THcAerogel("aero", "Aerogel");
+  SHMS->AddDetector(paero);
+  // Add calorimeter to SHMS apparatus
+  THcShower* pcal = new THcShower("cal", "Calorimeter");
+  SHMS->AddDetector(pcal);
+
+  // Add rastered beam apparatus
+  THaApparatus* pbeam = new THcRasteredBeam("P.rb", "Rastered Beamline");
+  gHaApps->Add(pbeam);
+
+  auto trkEff = new hcana::TrackingEfficiency("P.trackEfficiency", "SHMS tracking Efficiency", "P.hod");
+  gHaPhysics->Add(trkEff);
+  // Add physics modules
+  // Calculate reaction point
+  THaReactionPoint* prp = new THaReactionPoint("P.react", "SHMS reaction point", "P", "P.rb");
+  gHaPhysics->Add(prp);
+  // Calculate extended target corrections
+  THcExtTarCor* pext = new THcExtTarCor("P.extcor", "HMS extended target corrections", "P", "P.react");
+  gHaPhysics->Add(pext);
+  // Calculate golden track quantites
+  THaGoldenTrack* pgtr = new THaGoldenTrack("P.gtr", "SHMS Golden Track", "P");
+  gHaPhysics->Add(pgtr);
+  // Calculate the hodoscope efficiencies
+  THcHodoEff* peff = new THcHodoEff("phodeff", "SHMS hodo efficiency", "P.hod");
+  gHaPhysics->Add(peff);   
+
+  // Add event handler for scaler events
+  THcScalerEvtHandler* pscaler = new THcScalerEvtHandler("P", "Hall C scaler event type 1");
+  pscaler->AddEvtType(1);
+  pscaler->AddEvtType(4);
+  pscaler->AddEvtType(5);
+  pscaler->AddEvtType(6);
+  pscaler->AddEvtType(7);
+  pscaler->AddEvtType(129);
+  pscaler->SetDelayedType(129);
+  pscaler->SetUseFirstEvent(kTRUE);
+  gHaEvtHandlers->Add(pscaler);
+
+
+  // =================================
+  //  Kinematics Modules
+  // =================================
+
+  // Add Physics Module to calculate primary (scattered electrons) beam kinematics
+  //THcPrimaryKine* hkin_primary = new THcPrimaryKine("H.kin.primary", "HMS Single Arm Kinematics", "H", "H.rb");
+  //gHaPhysics->Add(hkin_primary);
+  // Add Physics Module to calculate secondary (scattered hadrons) beam kinematics
+  //THcSecondaryKine* pkin_secondary = new THcSecondaryKine("P.kin.secondary", "SHMS Single Arm Kinematics", "P", "H.kin.primary");
+  //gHaPhysics->Add(pkin_secondary);
+  THcPrimaryKine* kin = new THcPrimaryKine("P.kin", "SHMS Single Arm Kinematics", "P", "P.rb");
+  gHaPhysics->Add(kin);
+  
+  // =================================
+  //  Global Objects & Event Handlers
+  // =================================
+
+  // Add trigger apparatus
+  THaApparatus* TRG = new THcTrigApp("T", "TRG");
+  gHaApps->Add(TRG);
+  // Add trigger detector to trigger apparatus
+  //THcTrigDet* coin = new THcTrigDet("coin", "Coincidence Trigger Information");
+  // Suppress missing reference time warnings for these event types
+  //coin->SetEvtType(1);
+  //coin->AddEvtType(2);
+  //TRG->AddDetector(coin); 
+
+  // Add helicity detector to grigger apparatus
+  THcHelicity* helicity = new THcHelicity("helicity","Helicity Detector");
+  TRG->AddDetector(helicity);
+
+  // Add event handler for prestart event 125.
+  THcConfigEvtHandler* ev125 = new THcConfigEvtHandler("HC", "Config Event type 125");
+  gHaEvtHandlers->Add(ev125);
+  // Add event handler for EPICS events
+  THaEpicsEvtHandler* hcepics = new THaEpicsEvtHandler("epics", "HC EPICS event type 180");
+  gHaEvtHandlers->Add(hcepics);
+
+  // Add coin physics module 
+  //THcCoinTime* coinTime = new THcCoinTime("CTime", "Coincidende Time Determination", "P", "H", "T.coin");
+  //gHaPhysics->Add(coinTime);
+ 
+  // -----------------------------------------------------------
+  // Scandalizer 
+  // -----------------------------------------------------------
+  //
+  hcana::Scandalizer* analyzer = new hcana::Scandalizer;
+  //analyzer->_skip_events = 100;
+  analyzer->SetCodaVersion(2);
+  //analyzer->EnableBenchmarks(true);
+
+  // The following analyzes the first 2000 events (for pedestals, is required) 
+  // then  repeatedly skips 3000 events and processes 1000.
+  //auto pp0 = new hallc::scandalizer::SkipPeriodicAfterPedestal();
+  auto pp0 = new hallc::scandalizer::SkipAfterPedestal();
+  pp0->_analyzer = analyzer;
+
+  //SimplePostProcess([&]() { return 0; },
+  //                                                     [&](const THaEvData* evt) {
+  //                                                       static int counter = 0;
+  //                                                       if (evt->GetEvNum() > 2000) {
+  //                                                         if (counter == 0) {
+  //                                                           analyzer->_skip_events = 3000;
+  //                                                           counter                = 1000;
+  //                                                         } else {
+  //                                                           counter--;
+  //                                                         }
+  //                                                       }
+  //                                                       return 0;
+  //                                                     });
+
+  //hallc::scandalizer::SpectrometerMonitor * pp1a = new hallc::scandalizer::SpectrometerMonitor(hhod,hcer,hdc);
+  //pp1a->_analyzer = analyzer;
+  //pp1a->_spectrometer_name = "HMS";
+
+  hallc::scandalizer::SpectrometerMonitor * pp1 = new hallc::scandalizer::SpectrometerMonitor(phod,phgcer,pdc);
+  pp1->_analyzer = analyzer;
+
+  //hallc::scandalizer::SimplePostProcess* pp1 = new hallc::scandalizer::SimplePostProcess(
+  //  [&](){
+  //    return 0;
+  //  },
+  //  [&](const THaEvData* evt){
+  //    static int counter = 0;
+  //    static double eff_num             = 0.0000001;
+  //    static double eff_den             = 0.0;
+  //    static int    n_num               = 0;
+  //    static int    n_den               = 0;
+  //    int           shmsDC1Planes_nhits = 0;
+  //    int           shmsDC2Planes_nhits = 0;
+  //    for (int ip = 0; ip < 6; ip++) {
+  //      shmsDC1Planes_nhits += pdc->GetPlane(ip)->GetNHits();
+  //    }
+  //    for (int ip = 6; ip < 12; ip++) {
+  //      shmsDC2Planes_nhits += pdc->GetPlane(ip)->GetNHits();
+  //    }
+  //    bool   shms_DC_too_many_hits = (shmsDC1Planes_nhits > 8) || (shmsDC2Planes_nhits > 8);
+  //    double beta                  = phod->GetBetaNotrk();
+  //    bool   good_beta             = beta > 0.4;
+  //    bool   shms_good_hodoscope   = phod->fGoodScinHits;
+  //    bool   good_ntracks          = (pdc->GetNTracks() > 0);
+  //    bool   good_hgc              = phgcer->GetCerNPE() > 1;
+  //    if ((good_beta && shms_good_hodoscope) && (!shms_DC_too_many_hits) && good_hgc) {
+  //      eff_den = eff_den + 1.0;
+  //      n_den++;
+  //      if (good_ntracks) {
+  //        eff_num = eff_num + 1.0;
+  //        n_num++;
+  //      }
+  //    }
+  //    if ((evt->GetEvNum() > 1200) && (counter > 1000)) {
+  //      std::cout << " efficiency :  " << eff_num / eff_den << "\n";
+  //      //std::cout << " Event : " << evt->GetEvNum() << "  ( " << evt->GetEvType() << ")\n";
+  //      pv_list.Put("hcSHMSTrackingEff", eff_num/eff_den);
+  //      pv_list.Put("hcSHMSTrackingEff:Unc",std::sqrt(double(n_num))/(n_num+n_den+0.0000001));
+  //      pv_list.Put("hcSHMSTrackingEff.LOW",0.95);
+  //      pv_list.Put("hcSHMSTrackingEff.LOLO",0.93);
+
+  //      eff_num                = 0.000000001;
+  //      eff_den                = 0.0;
+  //      //analyzer->_skip_events = 300;
+  //      counter = 0;
+  //    }
+  //    counter++;
+  //    return 0; 
+  //  });
+
+  analyzer->AddPostProcess(pp0);
+  analyzer->AddPostProcess(pp1);
+  //analyzer->AddPostProcess(pp1a);
+
+  // A simple event class to be output to the resulting tree.
+  // Creating your own descendant of THaEvent is one way of
+  // defining and controlling the output.
+  THaEvent* event = new THaEvent;
+
+  // Define the run(s) that we want to analyze.
+  // We just set up one, but this could be many.
+  THcRun* run = new THcRun( pathList, Form(RunFileNamePattern, RunNumber) );
+
+  // Set to read in Hall C run database parameters
+  run->SetRunParamClass("THcRunParameters");
+  
+  // Eventually need to learn to skip over, or properly analyze the pedestal events
+  run->SetEventRange(1, MaxEvent); // Physics Event number, does not include scaler or control events.
+  run->SetNscan(1);
+  run->SetDataRequired(0x7);
+  //run->Print();
+
+  // Define the analysis parameters
+  TString ROOTFileName = Form(ROOTFileNamePattern, RunNumber, MaxEvent);
+  analyzer->SetCountMode(2);  // 0 = counter is # of physics triggers
+                              // 1 = counter is # of all decode reads
+                              // 2 = counter is event number
+
+  analyzer->SetEvent(event);
+  // Set EPICS event type
+  analyzer->SetEpicsEvtType(180);
+  // Define crate map
+  analyzer->SetCrateMapFileName("MAPS/db_cratemap.dat");
+  // Define output ROOT file
+  analyzer->SetOutFile(ROOTFileName.Data());
+  // Define DEF-file+
+  analyzer->SetOdefFile("DEF-files/SHMS/PRODUCTION/pstackana_production_all.def");
+  analyzer->SetOdefFile("DEF-files/SHMS/PRODUCTION/pstackana_production.def");
+  // Define cuts file
+  analyzer->SetCutFile("DEF-files/SHMS/PRODUCTION/CUTS/pstackana_production_cuts.def");  // optional
+  // File to record accounting information for cuts
+  analyzer->SetSummaryFile(Form("REPORT_OUTPUT/SHMS/PRODUCTION/summary_all_production_%d_%d.report", RunNumber, MaxEvent));  // optional
+  // Start the actual analysis.
+  analyzer->Process(run);
+  // Create report file from template
+  analyzer->PrintReport("TEMPLATES/SHMS/PRODUCTION/pstackana_production.template",
+  			Form("REPORT_OUTPUT/SHMS/PRODUCTION/replay_shms_all_production_%d_%d.report", RunNumber, MaxEvent));  // optional
+
+
+  delete analyzer;
+}