diff --git a/rich/rich_reco.py b/rich/options/rich_reco.py
similarity index 100%
rename from rich/rich_reco.py
rename to rich/options/rich_reco.py
diff --git a/rich/rich_hadrons.sh b/rich/rich_hadrons.sh
new file mode 100644
index 0000000000000000000000000000000000000000..c257059e54ba09c6e18724786d4d89d69080fb97
--- /dev/null
+++ b/rich/rich_hadrons.sh
@@ -0,0 +1,57 @@
+#!/bin/bash
+
+if [[ ! -n  "${JUGGLER_DETECTOR}" ]] ; then
+  export JUGGLER_DETECTOR="topside"
+fi
+
+if [[ ! -n  "${JUGGLER_N_EVENTS}" ]] ; then
+  export JUGGLER_N_EVENTS=100
+fi
+
+if [[ ! -n  "${E_start}" ]] ; then
+  export E_start=0.0
+fi
+
+if [[ ! -n  "${E_end}" ]] ; then
+  export E_end=30.0
+fi
+
+export JUGGLER_FILE_NAME_TAG="emcal_uniform_electrons"
+export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
+
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
+
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
+
+
+# Build the detector constructors.
+git clone https://eicweb.phy.anl.gov/EIC/detectors/${JUGGLER_DETECTOR}.git
+mkdir ${JUGGLER_DETECTOR}/build
+pushd ${JUGGLER_DETECTOR}/build
+cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j30 install
+popd
+
+# generate the input events
+# note datasets is now only used to develop datasets.
+#git clone https://eicweb.phy.anl.gov/EIC/datasets.git datasets
+root -b -q "ecal/scripts/emcal_electrons.cxx(${JUGGLER_N_EVENTS}, ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+root -b -q "ecal/scripts/emcal_electrons_reader.cxx(${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+
+pushd ${JUGGLER_DETECTOR}
+ls -l
+# run geant4 simulations
+python options/ForwardRICH/simu.py \
+       --compact=${JUGGLER_DETECTOR}.xml \
+       -i ../${JUGGLER_FILE_NAME_TAG} \
+       -o ${JUGGLER_SIM_FILE}
+       -n ${JUGGLER_N_EVENTS}
+
+# Need to figure out how to pass file name to juggler from the commandline
+xenv -x /usr/local/Juggler.xenv gaudirun.py ../ecal/options/rich_reco.py
+ls -l
+popd
+
+#@TODO add analysis scripts
+
diff --git a/rich/scripts/rich_samples.py b/rich/scripts/rich_samples.py
new file mode 100644
index 0000000000000000000000000000000000000000..0dcc2bdbcd5e73e31bb70fb4e0028dc37cf1112f
--- /dev/null
+++ b/rich/scripts/rich_samples.py
@@ -0,0 +1,44 @@
+import ROOT
+import numpy as np
+from ROOT import gROOT
+import argparse
+from matplotlib import pyplot as plt
+from matplotlib.patches import Ellipse
+import matplotlib.transforms as transforms
+
+parser = argparse.ArgumentParser(description='rich reconstruction analysis')
+parser.add_argument('nevents', type=int, help='number of events')
+args = parser.parse_args()
+
+gROOT.Macro('rootlogon.C')
+f = ROOT.TFile("rich_test_reco.root")
+events = f.events
+
+nev = 0
+while nev < args.nevents:
+    iev = np.random.randint(0, events.GetEntries())
+    events.GetEntry(iev)
+    if events.RICHClusters.size() == 0:
+        continue
+
+    nev += 1
+    x, y = [], []
+    for hit in events.RecoForwardRICHHits:
+        x.append(hit.position.x)
+        y.append(hit.position.y)
+
+    fig, ax = plt.subplots(figsize=(9, 9), dpi=120)
+    colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
+    ax.scatter(x, y, s=9, marker='s', color=colors[0])
+
+    for cluster in events.RICHClusters:
+        ellipse = Ellipse((cluster.position.x, cluster.position.y), width=cluster.radius * 2, height=cluster.radius * 2,
+                          edgecolor='red', fill=False, lw=2, alpha=0.8)
+        ellipse.set_transform(ax.transData)
+        ax.add_patch(ellipse)
+    ax.set_xlim(-60, 60)
+    ax.set_ylim(-60, 60)
+    ax.set_xlabel('X (mm)')
+    ax.set_ylabel('Y (mm)')
+    fig.savefig('rich_cluster_{}.png'.format(iev))
+