diff --git a/benchmarks/LGC/optical_photon_count.py b/benchmarks/LGC/optical_photon_count.py
new file mode 100644
index 0000000000000000000000000000000000000000..1916120960461d9c6f3542b7c96fdd7a84709eb5
--- /dev/null
+++ b/benchmarks/LGC/optical_photon_count.py
@@ -0,0 +1,68 @@
+import os
+import argparse
+import awkward as ak
+import uproot
+import numpy as np
+from matplotlib import pyplot as plt
+
+
+mc_cols = [
+    'MCParticles.generatorStatus',
+    'MCParticles.momentum.x',
+    'MCParticles.momentum.y',
+    'MCParticles.momentum.z',
+    'MCParticles.PDG'
+]
+
+lgc_cols = [
+    'LightGasCherenkovHits.cellID',
+    'LightGasCherenkovHits.time',
+]
+
+qe = 0.15
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description='Count optical photons on the PMT Arrays')
+    parser.add_argument('sim_file', type=str,
+                        help='path to simulation output (root file).')
+    parser.add_argument('-o', '--output-dir', type=str,
+                        default='.',
+                        help='output directory.')
+    parser.add_argument('-n', dest='nevents', type=int,
+                        default=-1,
+                        help='number of event to process.')
+    parser.add_argument('-s', '--entry-start', type=int,
+                        default=-1,
+                        help='event number to start')
+    args = parser.parse_args()
+
+    # determine event range
+    entry_start = None
+    entry_stop = None
+    if args.nevents > 0:
+        entry_stop = args.nevents
+    if args.entry_start > 0:
+        entry_start = args.entry_start
+        entry_stop += entry_start
+
+    events = uproot.open(args.sim_file)['events']
+
+    # generated particles
+    mc_pars = events.arrays(mc_cols, entry_stop=entry_stop, library='ak')
+    first_par_mask = mc_pars['MCParticles.generatorStatus'] == 1
+    mc_pars = mc_pars[first_par_mask]
+
+    momentum = np.sqrt(mc_pars['MCParticles.momentum.x']**2 + mc_pars['MCParticles.momentum.y']**2 + mc_pars['MCParticles.momentum.z']**2)
+    phi = ak.flatten(np.arctan2(mc_pars['MCParticles.momentum.x'], mc_pars['MCParticles.momentum.y']))
+    theta = ak.flatten(np.arccos(mc_pars['MCParticles.momentum.z']/momentum))
+
+    # LGC readouts
+    lgc_hits = events.arrays(lgc_cols, entry_stop=entry_stop, library='ak')
+    raw_counts = ak.count(lgc_hits['LightGasCherenkovHits.cellID'], axis=-1)
+    qe_counts = raw_counts*qe
+    counts_err = np.sqrt(qe * (1. - qe) * raw_counts)
+
+    # plots
+    fig, ax = plt.subplots(figsize=(8, 6), dpi=160)
+    ax.hist(qe_counts, bins=np.linspace(0, 60, 61), ec='k')
+    fig.savefig(os.path.join(args.output_dir, 'oph_counts.png'))
diff --git a/benchmarks/LGC/steering.py b/benchmarks/LGC/steering.py
index ffb39ed65759a7833af43276c1290606fc864750..37f67c743e80eaf19d4983cef38dd2055fe15142 100644
--- a/benchmarks/LGC/steering.py
+++ b/benchmarks/LGC/steering.py
@@ -1,5 +1,6 @@
 """
 A steering file that enables optical photon counting for LGC
+It sets up physics lists and particle filters for LGC, as well as setting a particle gun for SoLID configuration
    @author  C.Peng
 
 [simulation]: