diff --git a/benchmarks/far_forward/analysis/analysis_zdc_sipmontile_particles.cxx b/benchmarks/far_forward/analysis/analysis_zdc_sipmontile_particles.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..e135834514452179578265c1c4e80b4aacb6eec8
--- /dev/null
+++ b/benchmarks/far_forward/analysis/analysis_zdc_sipmontile_particles.cxx
@@ -0,0 +1,307 @@
+////////////////////////////////////////
+// Read ROOT output file
+// Plot variables
+// Jihee Kim 09/2021
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+#include <vector>
+
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/SimCalorimeterHitCollection.h"
+#include "edm4eic/ClusterCollection.h"
+#include "edm4eic/ClusterData.h"
+
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TGaxis.h"
+#include "TPaveStats.h"
+#include "TMath.h"
+#include "TH1.h"
+#include "TF1.h"
+#include "TH1D.h"
+#include "TVector3.h"
+
+R__LOAD_LIBRARY(libedm4eic.so)
+#include "fmt/core.h"
+#include "DD4hep/Detector.h"
+#include "DDG4/Geant4Data.h"
+#include "DDRec/CellIDPositionConverter.h"
+
+using ROOT::RDataFrame;
+
+void analysis_zdc_sipmontile_particles(
+  const std::string& input_fname = "sim_zdc_sipmontile_uniform.edm4hep.root",
+  const std::string& results_path = "results/far_forward/zdc_sipmontile/"
+) {
+  // Setting for graphs
+  TStyle* kStyle = new TStyle("kStyle","Kim's Style");
+  kStyle->SetOptStat("emr");
+  kStyle->SetOptTitle(0);
+  kStyle->SetOptFit(1101);
+  kStyle->SetStatColor(0);
+  kStyle->SetStatW(0.15);
+  kStyle->SetStatH(0.10);
+  kStyle->SetStatX(0.9);
+  kStyle->SetStatY(0.9);
+  kStyle->SetStatBorderSize(1);
+  kStyle->SetLabelSize(0.045,"xyz");
+  kStyle->SetTitleSize(0.045,"xyz");
+  kStyle->SetTitleOffset(1.2,"y");
+  kStyle->SetLineWidth(2);
+  kStyle->SetTitleFont(42,"xyz");
+  kStyle->SetLabelFont(42,"xyz");
+  kStyle->SetCanvasDefW(700);
+  kStyle->SetCanvasDefH(500);
+  kStyle->SetCanvasColor(0);
+  kStyle->SetPadTickX(1);
+  kStyle->SetPadTickY(1);
+  kStyle->SetPadGridX(1);
+  kStyle->SetPadGridY(1);
+  kStyle->SetPadLeftMargin(0.1);
+  kStyle->SetPadRightMargin(0.1);
+  kStyle->SetPadTopMargin(0.1);
+  kStyle->SetPadBottomMargin(0.1);
+  TGaxis::SetMaxDigits(3);
+  gROOT->SetStyle("kStyle");
+
+  ROOT::EnableImplicitMT();
+  ROOT::RDataFrame d0("events", input_fname);
+
+  // Thrown Energy [GeV]
+  auto E_gen = [](const std::vector<edm4hep::MCParticleData> & input) {
+    auto p = input[2];
+    auto energy = TMath::Sqrt(p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y + p.momentum.z * p.momentum.z + p.mass * p.mass);
+    return energy;
+  };
+
+  // Theta [rad]
+  auto Theta_gen = [](const std::vector<edm4hep::MCParticleData> & input) {
+    auto p = input[2];
+    auto theta = std::atan2(std::hypot(p.momentum.x, p.momentum.y), p.momentum.z);
+    return theta;
+  };
+  auto Thetarec = [](const std::vector<edm4eic::ClusterData> & input) {
+    std::vector<float> output;
+    for (const auto &c: input) {
+      auto r = c.position;
+      auto theta = std::atan2(std::hypot(r.x, r.y), r.z);
+      output.push_back(theta);
+    }
+    return output;
+  };
+
+  // Phi [rad]
+  auto Phi_gen = [](const std::vector<edm4hep::MCParticleData> & input) {
+    auto p = input[2];
+    auto phi = std::atan2(p.momentum.y, p.momentum.x);
+    if (phi < 0) phi += 2.*M_PI;
+    return phi;
+  };
+  auto Phirec = [](const std::vector<edm4eic::ClusterData> & input) {
+    std::vector<float> output;
+    for (const auto &c: input) {
+      auto r = c.position;
+      auto phi = std::atan2(r.y, r.x);
+      if (phi < 0) phi += 2.*M_PI;
+      output.push_back(phi);
+    }
+    return output;
+  };
+  
+  // Eta
+  auto Eta_gen = [](const std::vector<edm4hep::MCParticleData> & input) {
+    auto p = input[2];
+    auto theta = std::atan2(std::hypot(p.momentum.x, p.momentum.y), p.momentum.z);
+    auto eta = -std::log(std::tan(theta/2));
+    return eta;
+  };
+  auto Etarec = [](const std::vector<edm4eic::ClusterData> & input) {
+    std::vector<float> output;
+    for (const auto &c: input) {
+      auto r = c.position;
+      auto theta = std::atan2(std::hypot(r.x, r.y), r.z);
+      auto eta = -std::log(std::tan(theta/2));
+      output.push_back(eta);
+    }
+    return output;
+  };
+
+  // Define variables
+  auto d1 = d0
+    // Truth
+    .Define("E_gen",       E_gen,     {"MCParticles"})
+    .Define("Theta_gen",   Theta_gen, {"MCParticles"})
+    .Define("Phi_gen",     Phi_gen,   {"MCParticles"})
+    .Define("Eta_gen",     Eta_gen,   {"MCParticles"})
+    // Ecal
+    /*.Define("E_Ecal",      {"ZDCEcalClusters.energy"})
+    .Define("Theta_Ecal",  Thetarec, {"ZDCEcalClusters"})
+    .Define("Phi_Ecal",    Phirec,   {"ZDCEcalClusters"})
+    .Define("Eta_Ecal",    Etarec,   {"ZDCEcalClusters"})*/
+    // HCal
+    /*.Define("E_Hcal",      {"HcalFarForwardZDCClusters.energy"})
+    .Define("Theta_Hcal",  Thetarec, {"HcalFarForwardZDCClusters"})
+    .Define("Phi_Hcal",    Phirec,   {"HcalFarForwardZDCClusters"})
+    .Define("Eta_Hcal",    Etarec,   {"HcalFarForwardZDCClusters"})*/
+	;
+
+  // Histogram boundaries
+  const auto E_nbins = 100;
+  const auto E_max = 50.0;
+  const auto theta_nbins = 50;
+  const auto theta_max = 0.050;
+  const auto phi_nbins = 90;
+  const auto phi_min = 0.0;
+  const auto phi_max = 2.*M_PI;
+  const auto eta_nbins = 40;
+  const auto eta_min = 6.0;
+  const auto eta_max = 10.0;
+
+  // Define Histograms
+  auto hE_gen        = d1.Histo1D({"hE_gen",        "Thrown Energy; E_{gen} [GeV]; Events",                E_nbins, 0.0, E_max},  "E_gen");
+  auto hTheta_gen    = d1.Histo1D({"hTheta_gen",    "Thrown Theta; #theta_{gen} [rad]; Events",           theta_nbins, 0.0, theta_max},  "Theta_gen");
+  auto hPhi_gen      = d1.Histo1D({"hPhi_gen",      "Thrown Phi; #phi_{gen} [rad]; Events",                phi_nbins, phi_min, phi_max},  "Phi_gen");
+  auto hEta_gen      = d1.Histo1D({"hEta_gen",      "Thrown Eta; #eta_{gen}; Events",                      eta_nbins, eta_min, eta_max}, "Eta_gen");
+  auto hPhiTheta_gen = d1.Histo2D({"hPhiTheta_gen", "Thrown Phi vs Theta; #theta_{gen} [rad]; #phi_{gen}", theta_nbins, 0.0, theta_max, phi_nbins, phi_min, phi_max}, 
+                  "Theta_gen", "Phi_gen");
+  /*  auto hE_Ecal        = d1.Histo1D({"hE_Ecal",        "Reconstructed Ecal Energy; E_{rec} [GeV]; Events",                E_nbins, 0.0, E_max},  "E_Ecal");
+  auto hTheta_Ecal    = d1.Histo1D({"hTheta_Ecal",    "Reconstructed Ecal Theta; #theta_{rec} [rad]; Events",           theta_nbins, 0.0, theta_max},  "Theta_Ecal");
+  auto hPhi_Ecal      = d1.Histo1D({"hPhi_Ecal",      "Reconstructed Ecal Phi; #phi_{rec} [rad]; Events",                phi_nbins, phi_min, phi_max},  "Phi_Ecal");
+  auto hEta_Ecal      = d1.Histo1D({"hEta_Ecal",      "Reconstructed Ecal Eta; #eta_{rec}; Events",                      eta_nbins, eta_min, eta_max}, "Eta_Ecal");
+  auto hPhiTheta_Ecal = d1.Histo2D({"hPhiTheta_Ecal", "Reconstructed Ecal Phi vs Theta; #theta_{rec} [rad]; #phi_{rec}", theta_nbins, 0.0, theta_max, phi_nbins, phi_min, phi_max}, 
+  "Theta_Ecal", "Phi_Ecal");*/
+  /*
+  auto hE_Hcal        = d1.Histo1D({"hE_Hcal",        "Reconstructed Hcal Energy; E_{rec} [GeV]; Events",                E_nbins, 0.0, E_max},  "E_Hcal");
+  auto hTheta_Hcal    = d1.Histo1D({"hTheta_Hcal",    "Reconstructed Hcal Theta; #theta_{rec} [rad]; Events",           theta_nbins, 0.0, theta_max},  "Theta_Hcal");
+  auto hPhi_Hcal      = d1.Histo1D({"hPhi_Hcal",      "Reconstructed Hcal Phi; #phi_{rec} [rad]; Events",                phi_nbins, phi_min, phi_max},  "Phi_Hcal");
+  auto hEta_Hcal      = d1.Histo1D({"hEta_Hcal",      "Reconstructed Hcal Eta; #eta_{rec}; Events",                      eta_nbins, eta_min, eta_max}, "Eta_Hcal");
+  auto hPhiTheta_Hcal = d1.Histo2D({"hPhiTheta_Hcal", "Reconstructed Hcal Phi vs Theta; #theta_{rec} [rad]; #phi_{rec}", theta_nbins, 0.0, theta_max, phi_nbins, phi_min, phi_max}, 
+                  "Theta_Hcal", "Phi_Hcal");
+  */
+  // Event Counts
+  auto nevents_thrown = d1.Count();
+  std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
+
+  // Draw Histograms
+  {
+    TCanvas* c1 = new TCanvas("c1", "c1");
+    auto h = hE_gen->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c1->SaveAs(TString(results_path + "/zdc_E_gen.png"));
+    c1->SaveAs(TString(results_path + "/zdc_E_gen.pdf"));
+  }
+  {
+    TCanvas* c2 = new TCanvas("c2", "c2");
+    auto h = hTheta_gen->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c2->SaveAs(TString(results_path + "/zdc_Theta_gen.png"));
+    c2->SaveAs(TString(results_path + "/zdc_Theta_gen.pdf"));
+  }
+  {
+    TCanvas* c3 = new TCanvas("c3", "c3");
+    auto h = hPhi_gen->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c3->SaveAs(TString(results_path + "/zdc_Phi_gen.png"));
+    c3->SaveAs(TString(results_path + "/zdc_Phi_gen.pdf"));
+  }
+  {
+    TCanvas* c4 = new TCanvas("c4", "c4");
+    auto h = hEta_gen->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c4->SaveAs(TString(results_path + "/zdc_Eta_gen.png"));
+    c4->SaveAs(TString(results_path + "/zdc_Eta_gen.pdf"));
+  }
+  {
+    TCanvas* c5 = new TCanvas("c5", "c5");
+    auto h = hPhiTheta_gen->DrawCopy("COLZ");
+    c5->SaveAs(TString(results_path + "/zdc_PhiTheta_gen.png"));
+    c5->SaveAs(TString(results_path + "/zdc_PhiTheta_gen.pdf"));
+  }
+  /*
+  // Draw Histograms Ecal
+  {
+    TCanvas* c1 = new TCanvas("c1", "c1");
+    auto h = hE_Ecal->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c1->SaveAs(TString(results_path + "/zdc_E_Ecal.png"));
+    c1->SaveAs(TString(results_path + "/zdc_E_Ecal.pdf"));
+  }
+  {
+    TCanvas* c2 = new TCanvas("c2", "c2");
+    auto h = hTheta_Ecal->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c2->SaveAs(TString(results_path + "/zdc_Theta_Ecal.png"));
+    c2->SaveAs(TString(results_path + "/zdc_Theta_Ecal.pdf"));
+  }
+  {
+    TCanvas* c3 = new TCanvas("c3", "c3");
+    auto h = hPhi_Ecal->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c3->SaveAs(TString(results_path + "/zdc_Phi_Ecal.png"));
+    c3->SaveAs(TString(results_path + "/zdc_Phi_Ecal.pdf"));
+  }
+  {
+    TCanvas* c4 = new TCanvas("c4", "c4");
+    auto h = hEta_Ecal->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c4->SaveAs(TString(results_path + "/zdc_Eta_Ecal.png"));
+    c4->SaveAs(TString(results_path + "/zdc_Eta_Ecal.pdf"));
+  }
+  {
+    TCanvas* c5 = new TCanvas("c5", "c5");
+    auto h = hPhiTheta_Ecal->DrawCopy("COLZ");
+    c5->SaveAs(TString(results_path + "/zdc_PhiTheta_Ecal.png"));
+    c5->SaveAs(TString(results_path + "/zdc_PhiTheta_Ecal.pdf"));
+  }
+  */
+  /*
+  // Draw Histograms Hcal
+  {
+    TCanvas* c1 = new TCanvas("c1", "c1");
+    auto h = hE_Hcal->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c1->SaveAs(TString(results_path + "/zdc_E_Hcal.png"));
+    c1->SaveAs(TString(results_path + "/zdc_E_Hcal.pdf"));
+  }
+  {
+    TCanvas* c2 = new TCanvas("c2", "c2");
+    auto h = hTheta_Hcal->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c2->SaveAs(TString(results_path + "/zdc_Theta_Hcal.png"));
+    c2->SaveAs(TString(results_path + "/zdc_Theta_Hcal.pdf"));
+  }
+  {
+    TCanvas* c3 = new TCanvas("c3", "c3");
+    auto h = hPhi_Hcal->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c3->SaveAs(TString(results_path + "/zdc_Phi_Hcal.png"));
+    c3->SaveAs(TString(results_path + "/zdc_Phi_Hcal.pdf"));
+  }
+  {
+    TCanvas* c4 = new TCanvas("c4", "c4");
+    auto h = hEta_Hcal->DrawCopy();
+    h->SetLineWidth(3);
+    h->SetLineColor(kBlue);
+    c4->SaveAs(TString(results_path + "/zdc_Eta_Hcal.png"));
+    c4->SaveAs(TString(results_path + "/zdc_Eta_Hcal.pdf"));
+  }
+  {
+    TCanvas* c5 = new TCanvas("c5", "c5");
+    auto h = hPhiTheta_Hcal->DrawCopy("COLZ");
+    c5->SaveAs(TString(results_path + "/zdc_PhiTheta_Hcal.png"));
+    c5->SaveAs(TString(results_path + "/zdc_PhiTheta_Hcal.pdf"));
+  }*/
+}
diff --git a/benchmarks/far_forward/analysis/gen_zdc_sipmontile_particles.cxx b/benchmarks/far_forward/analysis/gen_zdc_sipmontile_particles.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..314416337bb46a7cd69c39728d113b3a1406fd0d
--- /dev/null
+++ b/benchmarks/far_forward/analysis/gen_zdc_sipmontile_particles.cxx
@@ -0,0 +1,97 @@
+//////////////////////////////////////////////////////////////
+// ZDC detector
+// Single Particle dataset
+// J. Kim 09/2021
+//////////////////////////////////////////////////////////////
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/Print.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+
+#include <TMath.h>
+#include <cmath>
+#include <iostream>
+#include <math.h>
+#include <random>
+#include <TRandom.h>
+#include <Math/Vector3D.h>
+#include <Math/RotationY.h>
+
+using namespace HepMC3;
+
+void gen_zdc_sipmontile_particles(int n_events = 1e6, std::string particle = "neutron", double p_start = 20.0, double p_end = 140.0, const std::string& out_fname = "zdc_neutron.hepmc") {
+  WriterAscii hepmc_output(out_fname);
+  int events_parsed = 0;
+  GenEvent evt(Units::GEV, Units::MM);
+
+  // Random number generator
+  TRandom* r1 = new TRandom();
+
+  // Set crossing angle [rad]
+  double crossing_angle = -0.025;
+
+  // Set polar angle range [rad]
+  double theta_min = 0.0;
+  double theta_max = 0.005;
+
+  for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
+    // FourVector(px,py,pz,e,pdgid,status)
+    // status - type 4 is beam
+    // pdgid 11 - electron
+    // pdgid 2212 - proton
+    GenParticlePtr p1 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
+    GenParticlePtr p2 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
+
+    // Define momentum
+    Double_t p_mag    = r1->Uniform(p_start, p_end);
+    Double_t phi      = r1->Uniform(0.0, 2.0 * M_PI);
+    Double_t theta    = acos(r1->Uniform(cos(theta_max), cos(theta_min)));
+    Double_t p0_x     = p_mag * std::cos(phi) * std::sin(theta);
+    Double_t p0_y     = p_mag * std::sin(phi) * std::sin(theta);
+    Double_t p0_z     = p_mag * std::cos(theta);
+
+    // Rotate the vector in Y by crossing angle when particles are being generated
+    ROOT::Math::XYZVector p0{p0_x, p0_y, p0_z};
+    ROOT::Math::RotationY r(crossing_angle);
+    auto p = r * p0;
+    auto px = p.X();
+    auto py = p.Y();
+    auto pz = p.Z();
+
+    // FourVector(px,py,pz,e,pdgid,status)
+    // status - type 1 is final state
+    // pdgid 22 - photon massless
+    // pdgid 2112 - neutron 939.565 MeV/c^2
+    GenParticlePtr p3;
+    if (particle == "neutron") {
+      p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, std::hypot(p_mag, 0.939565)), 2112, 1);
+    } else if (particle == "photon") {
+      p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, p_mag), 22, 1);
+    } else if (particle == "muon") {
+      p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, std::hypot(p_mag, 0.1057)), 13, 1);
+    } else {
+      std::cout << "Particle type " << particle << " not recognized!" << std::endl;
+      exit(-1);
+    }
+
+    // Create vertex
+    GenVertexPtr v1 = std::make_shared<GenVertex>();
+    v1->add_particle_in(p1);
+    v1->add_particle_in(p2);
+    v1->add_particle_out(p3);
+    evt.add_vertex(v1);
+
+    if (events_parsed == 0) {
+      std::cout << "First event: " << std::endl;
+      Print::listing(evt);
+    }
+
+    hepmc_output.write_event(evt);
+    if (events_parsed % 100 == 0) {
+      std::cout << "Event: " << events_parsed << std::endl;
+    }
+    evt.clear();
+  }
+  hepmc_output.close();
+  std::cout << "Events parsed and written: " << events_parsed << std::endl;
+}
diff --git a/benchmarks/far_forward/analysis/make_plots_zdc_sipmontile_hexplit.py b/benchmarks/far_forward/analysis/make_plots_zdc_sipmontile_hexplit.py
new file mode 100644
index 0000000000000000000000000000000000000000..ce8a5e1b322586f76d84db7aba5afb165f4aa516
--- /dev/null
+++ b/benchmarks/far_forward/analysis/make_plots_zdc_sipmontile_hexplit.py
@@ -0,0 +1,29 @@
+import numpy as np, pandas as pd, matplotlib.pyplot as plt
+import sys
+
+infile=sys.argv[1]
+outdir=sys.argv[2]
+
+df=pd.read_csv(infile)
+print(df.columns)
+#first the px over pz
+plt.hist(df.eval("(x_reco_rw/z_reco_rw-x_truth/z_truth)*1000"), bins=np.linspace(-2,2, 20), histtype='step', label='HEXPLIT')
+plt.hist(df.eval("(x_reco/z_reco-x_truth/z_truth)*1000"), bins=np.linspace(-2,2, 20), histtype='step', label='baseline')
+plt.xlabel("$x_{\\rm reco}/z_{\\rm reco}-x_{\\rm truth}/z_{\\rm reco}$ [mrad]")
+plt.ylabel("events")
+plt.legend()
+print("saving to file,",outdir+"/residuals_x.png")
+plt.savefig(outdir+"/residuals_x.png")
+
+
+plt.figure()
+plt.hist(df.eval("(y_reco_rw/z_reco_rw-y_truth/z_truth)*1000"), bins=np.linspace(-2, 2, 20), histtype='step', label='HEXPLIT')
+plt.hist(df.eval("(y_reco/z_reco-y_truth/z_truth)*1000"), bins=np.linspace(-2, 2, 20), histtype='step', label='baseline')
+plt.xlabel("$y_{\\rm reco}/z_{\\rm reco}-y_{\\rm truth}/z_{\\rm reco}$ [mrad]")
+plt.ylabel("events")
+plt.legend()
+print("saving to file,",outdir+"/residuals_y.png")
+plt.savefig(outdir+"/residuals_y.png")
+
+print([np.std(df[a][abs(df[a])<70]) for a in "dx dx_rw dy dy_rw".split()])
+
diff --git a/benchmarks/far_forward/config.yml b/benchmarks/far_forward/config.yml
index 34008f9f09bb3f50ac61722544cb4ed617d9c7e8..f8c6e4b395a9ea10fdc9afb08cb692163f6cd11e 100644
--- a/benchmarks/far_forward/config.yml
+++ b/benchmarks/far_forward/config.yml
@@ -17,12 +17,27 @@ ZDC_far_forward:
   stage: run
   needs: ["far_forward:compile"]
   script:
+    - print_env.sh
     - bash benchmarks/far_forward/run_zdc.sh --particle $PARTICLE
   parallel:
     matrix:
       - PARTICLE: ["neutron", "photon"]
   allow_failure: true
 
+ZDC_sipmontile:
+  extends: .rec_benchmark
+  variables:
+    DETECTOR_CONFIG: epic_zdc_sipm_on_tile_only
+  stage: run
+  needs: ["far_forward:compile"]
+  script:
+    - print_env.sh
+    - bash benchmarks/far_forward/run_zdc_sipmontile.sh --particle $PARTICLE
+  parallel:
+    matrix:
+      - PARTICLE: ["neutron", "photon"]
+  allow_failure: true
+
 far_forward:collect:
   stage: collect
   needs:
diff --git a/benchmarks/far_forward/hexplit/hexplit.py b/benchmarks/far_forward/hexplit/hexplit.py
new file mode 100644
index 0000000000000000000000000000000000000000..ecbafe2b145cfebbe99cdd4ec10737c051c72ae6
--- /dev/null
+++ b/benchmarks/far_forward/hexplit/hexplit.py
@@ -0,0 +1,566 @@
+import uproot as ur, numpy as np, pandas as pd
+
+c_in_mm_per_ns = 299.792458
+
+#default values.  Can be configured
+Emin=0.000472*0.1 # 1/10 a MIP
+tmax=150 # 150 ns
+
+def get_xyzr_reco_no_reweighting(arrays, event, w0=4, weight_by_granularity=True, prefix="ZDC", MIP=0.00047):
+    x=arrays[f'{prefix}HitsReco.position.x'][event]
+    y=arrays[f'{prefix}HitsReco.position.y'][event]
+    z=arrays[f'{prefix}HitsReco.position.z'][event]
+    E=arrays[f'{prefix}HitsReco.energy'][event]
+    t=arrays[f'{prefix}HitsReco.time'][event] - z/c_in_mm_per_ns #correct for time of flight
+    sl=arrays[f'{prefix}HitsReco.dimension.x'][event]/2
+    #make selection of hits to be used
+    slc=(E>Emin) & (t<tmax)
+    x=x[slc]
+    y=y[slc]
+    z=z[slc]
+    E=E[slc]
+    sl=sl[slc]
+
+    if type(w0)!=float or w0 !=0 :
+        w=w0+np.log((E+.0000001)/sum(E))
+        w=w*(w>0)
+    else :
+        w=E
+    if weight_by_granularity:
+        w=w/sl**2
+    
+    sumw=np.sum(w+.0000001, axis=-1)
+    x_reco=np.sum(x*w, axis=-1)/sumw
+    y_reco=np.sum(y*w, axis=-1)/sumw
+    z_reco=np.sum(z*w, axis=-1)/sumw
+    r_reco=np.hypot(x_reco,y_reco)
+    return [x_reco,y_reco,z_reco,r_reco]
+
+def get_xyzr_truth(arrays, event, w0=4, weight_by_granularity=True, prefix="ZDC"):
+    #first determine the z value used in the recon
+    z=arrays[f'{prefix}HitsReco.position.z'][event]
+    E=arrays[f'{prefix}HitsReco.energy'][event]
+    t=arrays[f'{prefix}HitsReco.time'][event] - z/c_in_mm_per_ns #correct for time of flight
+    sl=arrays[f'{prefix}HitsReco.dimension.x'][event]
+    
+    #make selection
+    slc=(E>Emin) & (t<tmax)
+    z=z[slc]
+    E=E[slc]
+    t=t[slc]
+    sl=sl[slc]
+    
+    if type(w0)!= float or w0 !=0:
+        w=w0+np.log((E+.0000001)/sum(E))
+        w=w*(w>0)
+    else :
+        w=E
+    if weight_by_granularity:
+        w=w/sl**2
+    
+    z_reco=np.sum(z*w, axis=-1)/np.sum(w+.000000001, axis=-1)
+    
+    px=arrays["MCParticles.momentum.x"][event,2]
+    py=arrays["MCParticles.momentum.y"][event,2]
+    pz=arrays["MCParticles.momentum.z"][event,2]
+
+    x_truth=px/pz*z_reco
+    y_truth=py/pz*z_reco
+    z_truth=z_reco
+    r_truth=np.hypot(x_truth,y_truth)
+    
+    return [x_truth,y_truth,z_truth,r_truth]
+
+sqrt3=np.sqrt(3)
+sqrt5=np.sqrt(5)
+
+def mx(a,b):
+    return a*(a>b)+b*(b>=a)
+#with H3
+def get_xyzr_reco_reweighted_H3(arrays, event, w0=5, weight_by_granularity=True, prefix="ZDC", MIP=0.000472):
+    x=arrays[f'{prefix}HitsReco.position.x'][event]
+    y=arrays[f'{prefix}HitsReco.position.y'][event]
+    z=arrays[f'{prefix}HitsReco.position.z'][event] 
+    E=arrays[f'{prefix}HitsReco.energy'][event]
+    #lay=arrays['HcalEndcapPInsertHitsReco.layer'][event]
+    sl=arrays[f'{prefix}HitsReco.dimension.x'][event]/2
+    t=arrays[f'{prefix}HitsReco.time'][event]  - z/c_in_mm_per_ns #correct for time of flight
+    slc=(E>Emin) & (t<tmax)
+    x=x[slc]
+    y=y[slc]
+    z=z[slc]
+    E=E[slc]
+    sl=sl[slc]
+
+    minz=min(z)
+    dz=min(z[z!=minz])-minz
+
+    Etot=sum(E)
+    if type(w0) !=float or  w0 != 0.0:
+        thresh=Etot*np.exp(-np.max(w0))
+    else :
+        thresh = 0
+    #print("thresh", thresh)
+    
+    xnew=[]
+    ynew=[]
+    znew=[]
+    Enew=[]
+    slnew=[]
+    phi=np.linspace(0, np.pi*5/3, 6)
+    cph=np.cos(phi)
+    sph=np.sin(phi)
+    roll_cph=np.roll(cph,1)
+    roll_sph=np.roll(sph,1)
+    
+    #centers of the triangles minus that of the hexagon
+    dx_tri=(cph+roll_cph)/3
+    dy_tri=(sph+roll_sph)/3
+    
+    for i in range(len(x)):
+        if E[i]<thresh:
+            continue
+        neighbors_found=0
+        Eneighbors=[0,0,0,0,0,0]
+        for j in range(len(x)):
+
+            if abs(z[i]-z[j])>dz*1.1 or E[j]<thresh or z[i]==z[j]:
+                continue
+            dx=(x[j]-x[i])/sl[i]
+            dy=(y[j]-y[i])/sl[i]
+            if abs(dx)>1.1 or abs(dy)>1.1:
+                continue
+            tol=0.01
+
+            
+            if neighbors_found==6:
+                break
+            
+            for k in range(6):
+                if Eneighbors[k]:
+                    continue
+                if abs(dx-cph[k])<tol and abs(dy-sph[k])<tol:
+                    #print("found neighbor")
+                    Eneighbors[k]=E[j]
+                    #print(Eneighbors, E[j])
+                    neighbors_found+=1
+                    break
+        #a=thresh*b
+        
+        Eneighbors=np.array(Eneighbors)
+
+        reweight_energy=mx(Eneighbors,MIP)*mx(np.roll(Eneighbors,1),MIP) 
+        reweight_energy/=sum(reweight_energy)
+
+        for k in range(6):
+            if E[i]*reweight_energy[k]<thresh:
+                continue
+            xnew.append(x[i]+sl[i]*dx_tri[k])
+            ynew.append(y[i]+sl[i]*dy_tri[k])
+            znew.append(z[i])
+            Enew.append(E[i]*reweight_energy[k])
+            slnew.append(sl[i]/sqrt5)
+            
+    xnew=np.array(xnew)
+    ynew=np.array(ynew)
+    znew=np.array(znew)
+    Enew=np.array(Enew)
+    slnew=np.array(slnew)
+    
+    
+    if type(w0)!=float or  w0 !=0:
+        w=w0+np.log((Enew+.0000001)/Etot)
+        w=w*(w>0)
+    else :
+        w=Enew
+    if weight_by_granularity:
+        w=w/slnew**2
+
+    sumw=np.sum(w+.0000001, axis=-1)
+    x_reco=np.sum(xnew*w, axis=-1)/sumw
+    y_reco=np.sum(ynew*w, axis=-1)/sumw
+    z_reco=np.sum(znew*w, axis=-1)/sumw
+    r_reco=np.hypot(x_reco,y_reco)
+    return [x_reco,y_reco,z_reco,r_reco]
+
+#with H4                                                                                                     
+def get_xyzr_reco_reweighted_H4(arrays, event, w0=6, weight_by_granularity=True, prefix="ZDC", MIP=0.000472):
+    x=arrays[f'{prefix}HitsReco.position.x'][event]
+    y=arrays[f'{prefix}HitsReco.position.y'][event]
+    z=arrays[f'{prefix}HitsReco.position.z'][event]
+    E=arrays[f'{prefix}HitsReco.energy'][event]
+    t=arrays[f'{prefix}HitsReco.time'][event] - z/c_in_mm_per_ns #correct for time of flight
+    sl=arrays[f'{prefix}HitsReco.dimension.x'][event]/2
+    lx=arrays[f'{prefix}HitsReco.local.x'][event]
+    ly=arrays[f'{prefix}HitsReco.local.y'][event]
+    lz=arrays[f'{prefix}HitsReco.local.z'][event]
+    
+    # time and energy cuts
+    slc=(E>Emin) & (t<tmax)
+    x=x[slc]
+    y=y[slc]
+    z=z[slc]
+    E=E[slc]
+    sl=sl[slc]
+    
+    lx=lx[slc]
+    ly=ly[slc]
+    lz=lz[slc]
+
+    tol=0.01
+    #determine the spacing between layers
+    if len(lz)==0:
+        return None
+    minz=min(lz)
+    dz=min(lz[lz>minz+tol])-minz
+
+    Etot=sum(E)
+    if type(w0) !=float or  w0 != 0.0:
+        thresh=Etot*np.exp(-np.max(w0))
+    else :
+        thresh = 0
+
+
+    xnew=[]
+    ynew=[]
+    znew=[]
+    Enew=[]
+    slnew=[]
+    phi=np.linspace(0, np.pi*5/3, 6)
+    cph=np.cos(phi)
+    sph=np.sin(phi)
+
+    for i in range(len(x)):
+        if E[i]<thresh:
+            continue
+        # there are twelve neighboring cell positions where we need to determine the energy
+        Eneighbors=[0,0,0,0,0,0,0,0,0,0,0,0]
+        for j in range(len(x)):
+            if abs(lz[i]-lz[j])>dz*2.1 or  E[j]<thresh or  j == i :
+                continue
+            dx=(lx[j]-lx[i])/sl[i]
+            dy=(ly[j]-ly[i])/sl[i]
+            if abs(dx)>1.6 or abs(dy)>1.6:
+                continue
+            tol=0.01
+            for k in range(6):
+                if abs(dx-1.5*cph[k])<tol and abs(dy-1.5*sph[k])<tol:
+                    Eneighbors[k]+=E[j]
+                    break
+            for k in range(6):
+                if abs(dx+sqrt3/2*sph[k])<tol and abs(dy-sqrt3/2*cph[k])<tol:
+                    Eneighbors[k+6]+=E[j]
+                    break
+        #print(Eneighbors)
+        Eneighbors=mx(np.array(Eneighbors),MIP)
+        
+        
+        reweight_energy_1=Eneighbors[:6]*np.roll(Eneighbors[6:],1)*np.roll(Eneighbors[6:],2)
+        reweight_energy_2=Eneighbors[6:]*np.roll(Eneighbors[6:],-1)*np.roll(Eneighbors[6:],1)
+        reweight_energy=np.concatenate([reweight_energy_1, reweight_energy_2])
+        reweight_energy/=sum(reweight_energy)
+
+        for k in range(6):
+            if E[i]*reweight_energy[k]< thresh:
+                continue
+            xnew.append(x[i]+sl[i]*0.75*cph[k])
+            ynew.append(y[i]+sl[i]*0.75*sph[k])
+            znew.append(z[i])
+            Enew.append(E[i]*reweight_energy[k])
+            slnew.append(sl[i]/sqrt5)
+        
+        for k in range(6):
+            if E[i]*reweight_energy[k+6]<thresh:
+                continue
+            xnew.append(x[i]+sl[i]*sqrt3/4*-sph[k])
+            ynew.append(y[i]+sl[i]*sqrt3/4*cph[k])
+            znew.append(z[i])
+            Enew.append(E[i]*reweight_energy[k+6])
+            slnew.append(sl[i]/sqrt5)
+        
+    
+    xnew=np.array(xnew)
+    ynew=np.array(ynew)
+    znew=np.array(znew)
+    Enew=np.array(Enew)
+    
+    slnew=np.array(slnew)
+        
+    if type(w0)!=float or  w0 !=0:
+        w=w0+np.log((Enew+.0000001)/Etot)
+        w=w*(w>0)
+    else :
+        w=Enew
+    if weight_by_granularity:
+        w=w/slnew**2
+
+    sumw=np.sum(w+.0000001, axis=-1)
+    x_reco=np.sum(xnew*w, axis=-1)/sumw
+    y_reco=np.sum(ynew*w, axis=-1)/sumw
+    z_reco=np.sum(znew*w, axis=-1)/sumw
+    r_reco=np.hypot(x_reco,y_reco)
+
+    return [x_reco,y_reco,z_reco,r_reco]
+
+sqrt2=np.sqrt(2)
+#with S2
+def get_xyzr_reco_reweighted_S2(arrays, event, w0=6, weight_by_granularity=True, prefix="ZDC", MIP=0.000472):
+    x=arrays[f'{prefix}HitsReco.position.x'][event]
+    y=arrays[f'{prefix}HitsReco.position.y'][event]
+    z=arrays[f'{prefix}HitsReco.position.z'][event] 
+    E=arrays[f'{prefix}HitsReco.energy'][event]
+    
+    #side length.  
+    sl=arrays[f'{prefix}HitsReco.dimension.x'][event]
+    t=arrays[f'{prefix}HitsReco.time'][event] - z/c_in_mm_per_ns #correct for time of flight
+    slc=(E>Emin) & (t<tmax)
+    x=x[slc]
+    y=y[slc]
+    z=z[slc]
+    E=E[slc]
+    sl=sl[slc]
+
+    
+    minz=min(z)
+    dz=min(z[z!=minz])-minz
+
+    Etot=sum(E)
+    if type(w0) !=float or  w0 != 0.0:
+        thresh=Etot*np.exp(-np.max(w0))
+    else :
+        thresh = 0
+
+
+    xnew=[]
+    ynew=[]
+    znew=[]
+    Enew=[]
+    slnew=[]
+    phi=np.linspace(np.pi/4, np.pi*7/4, 4)
+    cph=np.cos(phi)
+    sph=np.sin(phi)
+    #print("bbb")
+    for i in range(len(x)):
+        if E[i]<thresh:
+            continue
+        neighbors_found = 0
+        Eneighbors=[0,0,0,0]
+        for j in range(len(x)):
+                                
+            if abs(z[i]-z[j])>dz*1.1 or  E[j]<thresh or  j == i or z[i]==z[j] :
+                continue
+            dx=(x[j]-x[i])/sl[i]
+            dy=(y[j]-y[i])/sl[i]
+            if abs(dx)>0.6 or abs(dy)>0.6:
+                continue
+            tol=0.01
+            for k in range(4):
+                if abs(dx-sqrt2*cph[k]/2)<tol and abs(dy-sqrt2*sph[k]/2)<tol:
+                    Eneighbors[k]+=E[j]
+                    neighbors_found+=1
+                    break
+            if neighbors_found == 8:
+                break
+            
+        #print("??")
+
+        Eneighbors=mx(np.array(Eneighbors),MIP)
+        #print("Eneighbors:", Eneighbors)
+        
+        reweight_energy=Eneighbors/sum(Eneighbors)
+        #print("reweight_energy:", reweight_energy)
+        #print(reweight_energy)
+        #print(Eneighbors)
+        #reweight_energy/=sum(reweight_energy)
+
+        #print("aaa")
+        for k in range(4):
+            if E[i]*reweight_energy[k]< thresh:
+                continue
+            xnew.append(x[i]+sl[i]*sqrt2/4*cph[k])
+            ynew.append(y[i]+sl[i]*sqrt2/4*sph[k])
+            znew.append(z[i])
+            Enew.append(E[i]*reweight_energy[k])
+            slnew.append(sl[i]/2)
+        
+    
+    xnew=np.array(xnew)
+    ynew=np.array(ynew)
+    znew=np.array(znew)
+    Enew=np.array(Enew)
+    
+    slnew=np.array(slnew)
+    #print("fff")
+        
+    if type(w0)!=float or  w0 !=0:
+        w=w0+np.log((Enew+.0000001)/Etot)
+        w=w*(w>0)
+    else :
+        w=Enew
+    if weight_by_granularity:
+        w=w/slnew**2
+    #print("ggg")
+    sumw=np.sum(w+.0000001, axis=-1)
+    x_reco=np.sum(xnew*w, axis=-1)/sumw
+    y_reco=np.sum(ynew*w, axis=-1)/sumw
+    z_reco=np.sum(znew*w, axis=-1)/sumw
+    r_reco=np.hypot(x_reco,y_reco)
+    #print("hh")
+    return [x_reco,y_reco,z_reco,r_reco]
+
+
+if __name__ == "__main__":
+    import argparse
+    parser = argparse.ArgumentParser(
+                    prog='process_showers',
+                    description='determines the position of the showers',
+                    epilog='Text at the bottom of help')
+    parser.add_argument('infile')
+    parser.add_argument('outfile')
+    parser.add_argument('--H4', help="flag for using H4 HEX-SPLIT",
+                    action='store_true')  # on/off flag
+    parser.add_argument('--H3',	help="flag for using H3 HEX-SPLIT",
+                    action='store_true')  # on/off flag 
+    parser.add_argument('--S2',  help="flag for using S2 HEX-SPLIT",
+                    action='store_true')  # on/off flag
+    parser.add_argument('-n', '--nevents', help="number of events to run", default=-1, type=int)
+    parser.add_argument("-p", '--prefix', help="prefix for the detector type (hit type)", default="ZDC")
+    parser.add_argument("-s", '--skip', help="number of events to skip", default=0, type=int)
+    parser.add_argument("--w0_rw", help="w0 used in reweighted", default=5, type=float)
+    parser.add_argument("--w0_nrw", help="w0 used in reweighted", default=4, type=float)
+    parser.add_argument("--w0_use_range", help="use a range for w0", action="store_true")
+    parser.add_argument("--MIP", help="MIP value in GeV", default=0.000472, type=float)
+    args = parser.parse_args()
+    
+    import sys
+    infile=args.infile
+    outfile=args.outfile
+    useH3reweighting=args.H3
+    useH4reweighting=args.H4
+    useS2reweighting=args.S2
+    nevents=args.nevents
+    prefix=args.prefix
+    first_event=args.skip
+    arrays=ur.open(f'{infile}:events').arrays()
+    w0=args.w0_rw
+    MIP=args.MIP
+    w0_nrw=args.w0_nrw
+
+    w0_use_range=args.w0_use_range
+    if w0_use_range:
+        w0_nrw=np.array([[a] for a in np.linspace(3.0, 8.0, 21)])
+        w0=np.array([[a] for a in np.linspace(3.0, 8.0, 21)])
+    r_truths=[]    
+    x_truths=[]
+    y_truths=[]
+    z_truths=[]
+
+    r_recos=[]
+    x_recos=[]
+    y_recos=[]
+    z_recos=[]
+
+    r_recos_rw=[]
+    x_recos_rw=[]
+    y_recos_rw=[]
+    z_recos_rw=[]
+    
+    drs=[]
+    drs_rw=[]
+    dxs=[]
+    dys=[]
+    dxs_rw=[]
+    dys_rw=[]
+    Es=[]
+    mc_pzs=[]
+    w0s=[]
+    if nevents==-1:
+        nevents=len(arrays)
+    
+    for event in range(first_event,first_event+nevents):
+        if len(arrays[f'{prefix}HitsReco.position.x'][event])== 0:
+            print("warning, empty event.  skipping")
+            continue
+        #print(arrays['ZDCHitsReco.energy'][event])
+        #print(len(arrays['ZDCHitsReco.position.x'][event]))
+        
+
+        x_reco, y_reco, z_reco, r_reco=get_xyzr_reco_no_reweighting(arrays, event, w0=w0_nrw, weight_by_granularity=True, prefix=prefix)
+        if useH3reweighting:
+            x_reco_rw, y_reco_rw, z_reco_rw, r_reco_rw=get_xyzr_reco_reweighted_H3(arrays, event, w0=w0, MIP=MIP, weight_by_granularity=True, prefix=prefix)
+        elif useH4reweighting:
+            ret=get_xyzr_reco_reweighted_H4(arrays, event, w0=w0, MIP=MIP, weight_by_granularity=True, prefix=prefix)
+            if ret is not None:
+                x_reco_rw, y_reco_rw, z_reco_rw, r_reco_rw=ret
+            else : continue
+        elif useS2reweighting:
+            x_reco_rw, y_reco_rw, z_reco_rw, r_reco_rw=get_xyzr_reco_reweighted_S2(arrays, event, w0=w0, MIP=MIP, weight_by_granularity=True, prefix=prefix)
+
+        x_truth, y_truth, z_truth, r_truth=get_xyzr_truth(arrays, event, w0=w0, weight_by_granularity=True, prefix=prefix)
+        #print(r_truth, r_reco)
+        drs.append(r_reco-r_truth)
+        dxs.append(x_reco-x_truth)
+        dys.append(y_reco-y_truth)
+        r_recos.append(r_reco)
+        x_recos.append(x_reco)
+        y_recos.append(y_reco)
+        z_recos.append(z_reco)
+        r_truths.append(r_truth)
+        x_truths.append(x_truth)
+        y_truths.append(y_truth)
+        z_truths.append(z_truth)
+        if useH3reweighting or useH4reweighting or useS2reweighting:
+            drs_rw.append(r_reco_rw-r_truth)
+            dxs_rw.append(x_reco_rw-x_truth)
+            dys_rw.append(y_reco_rw-y_truth)
+            r_recos_rw.append(r_reco_rw)
+            x_recos_rw.append(x_reco_rw)
+            y_recos_rw.append(y_reco_rw)
+            z_recos_rw.append(z_reco_rw)
+        Es.append(sum(arrays[f'{prefix}HitsReco.energy'][event]))
+        mc_pzs.append(arrays[f'MCParticles.momentum.z'][event,2])
+        if w0_use_range:
+            w0s.append(w0)
+        #except Exception as e:
+        #    print(e)
+        #     e#pass
+
+        if event%10==0:
+            print(f"{infile}: done with event {event}/{nevents}")
+    if not w0_use_range:
+        d=dict(E=Es, dr=drs, dy=dys, dx=dxs, mc_pz=mc_pzs, x_truth=x_truths, y_truth=y_truths, z_truth=z_truths,
+               x_reco=x_recos, y_reco=y_recos,z_reco=z_recos,r_reco=r_recos)
+        if useH3reweighting or useH4reweighting or useS2reweighting:
+            d["dr_rw"]=drs_rw
+            d["dx_rw"]=dxs_rw
+            d["dy_rw"]=dys_rw
+            d['r_reco_rw']=r_recos_rw
+            d['x_reco_rw']=x_recos_rw
+            d['y_reco_rw']=y_recos_rw
+            d['z_reco_rw']=z_recos_rw
+    else :
+        w0s = [a[0] for a in w0] # flatten array
+        d=dict(E=Es, dr=drs, dy=dys, dx=dxs, mc_pz=mc_pzs, x_truth=x_truths, y_truth=y_truths, z_truth=z_truths, r_truth=r_truths, x_reco=x_recos, y_reco=y_recos, z_reco=z_recos, r_reco=r_recos)
+        if useH3reweighting or useH4reweighting or useS2reweighting:
+            d["dr_rw"]=drs_rw
+            d["dx_rw"]=dxs_rw
+            d["dy_rw"]=dys_rw
+            d['r_reco_rw']=r_recos_rw
+            d['x_reco_rw']=x_recos_rw
+            d['y_reco_rw']=y_recos_rw
+            d['z_reco_rw']=z_recos_rw
+        for key in list(d.keys()):
+            #print(key)
+            #print(type(d[key]))
+            #print(d[key][0])
+            if key not in "E mc_pz w0s".split():
+                print("making new columns")
+                for i in range(len(w0s)):
+                    d[f"{key}_w0_{w0s[i]}".replace(".", "pt")]=[d[key][j][i] for j in range(len(d[key]))]
+            
+    
+    print({a:len(d[a]) for a in d})
+    if ".csv" in outfile:
+        pd.DataFrame(d).to_csv(outfile)
+    if ".pkl" in outfile:
+        pd.DataFrame(d).to_pickle(outfile)
diff --git a/benchmarks/far_forward/options/zdc_sipmontile_reconstruction.py b/benchmarks/far_forward/options/zdc_sipmontile_reconstruction.py
new file mode 100644
index 0000000000000000000000000000000000000000..bb986d452f4d621ea6c0a5fd92f9ca70d75ff8ce
--- /dev/null
+++ b/benchmarks/far_forward/options/zdc_sipmontile_reconstruction.py
@@ -0,0 +1,91 @@
+import os
+import ROOT
+from Gaudi.Configuration import *
+from GaudiKernel.SystemOfUnits import mm, MeV, rad, ns
+from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+
+from Configurables import PodioInput
+from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
+from Configurables import Jug__Reco__CalorimeterHitReco as CalHitReco
+from Configurables import Jug__Reco__CalorimeterHitsMerger as CalHitsMerge
+from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
+from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
+
+kwargs = dict()
+
+# input arguments through environment variables
+kwargs['input'] = os.environ.get('JUGGLER_SIM_FILE', '../epic/zdc_sipmontile_neutrons.edm4hep.root')
+kwargs['output'] = os.environ.get('JUGGLER_REC_FILE', 'rec_zdc_sipmontile_neutrons.root')
+kwargs['compact'] = os.environ.get('DETECTOR_COMPACT_PATH', '../epic/epic.xml')
+kwargs['nev'] = int(os.environ.get('JUGGLER_N_EVENTS', 100))
+kwargs['PbSci_sf'] = float(os.environ.get('ZDC_PbSCI_SAMP_FRAC', 1.0))
+
+if kwargs['nev'] < 1:
+    f = ROOT.TFile(kwargs['input'])
+    kwargs['nev'] = f.events.GetEntries()
+
+print(kwargs)
+
+geo_service = GeoSvc('GeoSvc', detectors=kwargs['compact'].split(','), OutputLevel=INFO)
+podioevent = EICDataSvc('EventDataSvc', inputs=kwargs['input'].split(','), OutputLevel=DEBUG)
+
+sim_colls = [
+    'MCParticles',
+    'HcalFarForwardZDCHits',
+    'HcalFarForwardZDCHitsContributions',
+
+]
+
+podin = PodioInput('PodioReader', collections=sim_colls, OutputLevel=DEBUG)
+podout = PodioOutput('out', filename=kwargs['output'])
+
+# Digitization
+
+ci_zdc_daq = dict(\
+         dynamicRangeADC=200.*MeV,
+         capacityADC=32768,
+         pedestalMean=400,
+         pedestalSigma=10)
+
+SiPMonTile_digi = CalHitDigi('SiPMonTile_digi',
+        inputHitCollection='HcalFarForwardZDCHits',
+                             outputHitCollection='HcalFarForwardZDCRawHits',
+                             **ci_zdc_daq)
+
+# Reconstruction
+SiPMonTile_reco = CalHitReco('SiPMonTile_reco',
+        inputHitCollection=SiPMonTile_digi.outputHitCollection,
+        outputHitCollection='HcalFarForwardZDCHitsReco',
+                             readoutClass='HcalFarForwardZDCHits',**ci_zdc_daq,
+                             thresholdFactor=0.0,
+                             samplingFraction="1.0")
+#        localDetFields=['system'])
+
+# Clustering
+'''SiPMonTile_cl = IslandCluster('SiPMonTile_cl',
+        inputHitCollection=SiPMonTile_reco.outputHitCollection,
+        outputProtoClusterCollection='HcalFarForwardZDCProtoClusters',
+        minClusterCenterEdep=3.*MeV,
+        minClusterHitEdep=0.1*MeV,
+        localDistXY=[50*mm, 50*mm],
+        splitCluster=True)
+
+SiPMonTile_clreco = RecoCoG('SiPMonTile_clreco',
+        inputProtoClusterCollection = SiPMonTile_cl.outputProtoClusterCollection,
+        outputClusterCollection='HcalFarForwardZDCClusters',
+        logWeightBase=3.6)
+'''
+podout.outputCommands = ['keep *']
+
+ApplicationMgr(
+    TopAlg=[podin,
+        SiPMonTile_digi, 
+        SiPMonTile_reco, 
+#        SiPMonTile_cl, SiPMonTile_clreco, 
+        podout],
+    EvtSel='NONE',
+    EvtMax=kwargs['nev'],
+    ExtSvc=[podioevent],
+    OutputLevel=DEBUG
+)
+
diff --git a/benchmarks/far_forward/run_zdc_sipmontile.sh b/benchmarks/far_forward/run_zdc_sipmontile.sh
new file mode 100644
index 0000000000000000000000000000000000000000..fc7d379c1e58f46fa705e7ba53609f5b5901b72c
--- /dev/null
+++ b/benchmarks/far_forward/run_zdc_sipmontile.sh
@@ -0,0 +1,154 @@
+#!/bin/bash
+
+function print_the_help {
+  echo "USAGE: ${0} [--rec-only] [--ana-only] [--particle=neutron] "
+  echo "  OPTIONS: "
+  echo "            --rec-only   only run gaudi reconstruction part"
+  echo "            --ana-only   only run analysis macros part"
+  exit 
+}
+
+PARTICLE="neutron"
+REC_ONLY=
+ANALYSIS_ONLY=
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+  key="$1"
+
+  case $key in
+    -h|--help)
+      shift # past argument
+      print_the_help
+      ;;
+    --rec-only)
+      REC_ONLY=1
+      shift # past value
+      ;;
+    --ana-only)
+      ANALYSIS_ONLY=1
+      shift # past value
+      ;;
+    --particle)
+      PARTICLE="$2"
+      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
+
+
+print_env.sh
+
+## To run the reconstruction, we need the following global variables:
+## - DETECTOR:         the detector package we want to use for this benchmark
+## - DETECTOR_VERSION: the detector package we want to use for this benchmark
+## - DETECTOR_PATH:            full path to the detector definitions
+##
+## You can ready options/env.sh for more in-depth explanations of the variables
+## and how they can be controlled.
+export DETECTOR_PATH=${DETECTOR_PATH}
+
+if [[ ! -n  "${JUGGLER_N_EVENTS}" ]] ; then 
+  export JUGGLER_N_EVENTS=100
+fi
+
+if [[ ! -n  "${E_start}" ]] ; then
+  export E_start=10.0
+fi
+
+if [[ ! -n  "${E_end}" ]] ; then
+  export E_end=10.0
+fi
+
+#if [[ ! -n "${ZDC_SCI_SAMP_FRAC}" ]] ; then
+#  export ZDC_PbSCI_SAMP_FRAC=1.0
+#fi
+
+export DETECTOR_CONFIG="epic_zdc_sipm_on_tile_only"
+echo "DETECTOR_CONFIG is set to "$DETECTOR_CONFIG
+
+export DETECTOR_COMPACT_PATH=${DETECTOR_PATH}/${DETECTOR_CONFIG}.xml
+
+export JUGGLER_FILE_NAME_TAG="zdc_sipmontile_${PARTICLE}"
+export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
+
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
+export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
+
+export RESULTS_PATH="results/far_forward/zdc_sipmontile/${PARTICLE}"
+export PHYSICS_LIST="FTFP_BERT"
+
+mkdir -p ${RESULTS_PATH}
+
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "DETECTOR = ${DETECTOR}"
+
+if [[ -z "${REC_ONLY}" && -z "${ANALYSIS_ONLY}" ]] ;
+then
+  echo "Generating input events"
+  root -b -q "benchmarks/far_forward/analysis/gen_zdc_sipmontile_particles.cxx(${JUGGLER_N_EVENTS}, \"${PARTICLE}\", ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running script: generating input events"
+    exit 1
+  fi
+
+  echo "Running Geant4 simulation"
+  npsim  \
+    --numberOfEvents ${JUGGLER_N_EVENTS} \
+    --physicsList ${PHYSICS_LIST} \
+    --compactFile ${DETECTOR_PATH}/${DETECTOR_CONFIG}.xml \
+    --inputFiles  ${JUGGLER_FILE_NAME_TAG}.hepmc \
+    --outputFile  ${JUGGLER_SIM_FILE}
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running npdet: Geant4 simulation"
+    exit 1
+  fi
+fi
+
+rootls -t ${JUGGLER_SIM_FILE}
+
+if [[ -z "${ANALYSIS_ONLY}" ]] ;
+then
+  gaudirun.py benchmarks/far_forward/options/zdc_sipmontile_reconstruction.py
+  status=$?
+  if [[ "$status" -ne "0" && "$status" -ne "4" ]] ; then
+    echo "ERROR running juggler, got $status"
+    exit 1
+  fi
+fi
+
+rootls -t ${JUGGLER_REC_FILE}
+
+HEXPLIT_OUT_FILE=${JUGGLER_REC_FILE}_hexplit.csv
+
+python benchmarks/far_forward/hexplit/hexplit.py ${JUGGLER_REC_FILE} ${HEXPLIT_OUT_FILE} --prefix HcalFarForwardZDC --w0_rw 5 --w0_nrw 5 --H4
+
+ls $HEXPLIT_OUT_FILE
+
+echo "Running analysis root scripts"
+root -b -q "benchmarks/far_forward/analysis/analysis_zdc_sipmontile_particles.cxx+(\"${JUGGLER_REC_FILE}\", \"${RESULTS_PATH}\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running root analysis script"
+  exit 1
+fi
+
+echo "making plots of hexplit residuals"
+python benchmarks/far_forward/analysis/make_plots_zdc_sipmontile_hexplit.py $HEXPLIT_OUT_FILE ${RESULTS_PATH}
+
+root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}")
+if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
+  # file must be less than 10 MB to upload
+  if [[ "${root_filesize}" -lt "10000000" ]] ; then
+    cp ${JUGGLER_SIM_FILE} results/far_forward/zdc_sipmontile/.
+    cp ${JUGGLER_REC_FILE} results/far_forward/zdc_sipmontile/.
+  fi
+fi