diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 5c314382ec1f7e2b339371ea6bd750c43dcc4452..a4ad26720e837d426627dbff932588bc21380cd9 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -88,7 +88,7 @@ include:
 
 final_report:
   stage: finish
-  needs: ["ecal_collect", "tracking_central_electrons", "clustering:results", "track_finding:collect", "track_fitting:collect"]
+  needs: ["ecal_collect", "tracking_central_electrons", "clustering:results", "track_finding:collect", "track_fitting:collect", "far_forward:collect"]
   script:
     # disabled while we address ACTS issues
     #- mkdir -p results/views && cd results/views && bash ../../bin/download_views
diff --git a/benchmarks/far_forward/analysis/analysis_zdc_neutrons.cxx b/benchmarks/far_forward/analysis/analysis_zdc_neutrons.cxx
deleted file mode 100644
index 810c443188c3762a3d265e0bcd7815de5e186e68..0000000000000000000000000000000000000000
--- a/benchmarks/far_forward/analysis/analysis_zdc_neutrons.cxx
+++ /dev/null
@@ -1,156 +0,0 @@
-////////////////////////////////////////
-// 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 "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(libeicd.so)
-#include "fmt/core.h"
-#include "DD4hep/Detector.h"
-#include "DDG4/Geant4Data.h"
-#include "DDRec/CellIDPositionConverter.h"
-
-using ROOT::RDataFrame;
-
-void analysis_zdc_neutrons(const char* input_fname = "sim_zdc_uniform_neutron.edm4hep.root")
-{
-  // 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 Ethr = [](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 [mrad]
-  auto Thetathr = [](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*1000.0;
-  };
-
-  // Phi [rad]
-  auto Phithr = [](const std::vector<edm4hep::MCParticleData> & input) {
-    auto p = input[2];
-    auto phi = std::atan2(p.momentum.y, p.momentum.x);
-    return phi;
-  };
-  
-  // Eta
-  auto Etathr = [](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;
-  };  
-
-  // Define variables
-  auto d1 = d0.Define("Ethr",      Ethr,     {"MCParticles"})
-              .Define("Thetathr",  Thetathr, {"MCParticles"})
-              .Define("Phithr",    Phithr,   {"MCParticles"})
-	      .Define("Etathr",    Etathr,   {"MCParticles"})
-	      ;
-
-  // Define Histograms
-  auto hEthr        = d1.Histo1D({"hEthr",        "Thrown Energy; E_{thr} [GeV]; Events",                100, 0.0, 0.3},  "Ethr");
-  auto hThetathr    = d1.Histo1D({"hThetathr",    "Thrown Theta; #theta_{thr} [mrad]; Events",           100, 0.0, 5.0},  "Thetathr");
-  auto hPhithr      = d1.Histo1D({"hPhithr",      "Thrown Phi; #phi_{thr} [rad]; Events",                100, 0.0, 5.0},  "Phithr");
-  auto hEtathr      = d1.Histo1D({"hEtathr",      "Thrown Eta; #eta_{thr}; Events",                      100, 6.0, 15.0}, "Etathr");
-  auto hPhiThetathr = d1.Histo2D({"hPhiThetathr", "Thrown Phi vs Theta; #theta_{thr} [mrad]; #phi_{thr}", 20, 0.0, 5.0, 20, 0.0, 5.0}, 
-                  "Thetathr", "Phithr");
-
-  // 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 = hEthr->DrawCopy();
-    h->SetLineWidth(3);
-    h->SetLineColor(kBlue);
-    c1->SaveAs("results/far_forward/zdc/zdc_neutron_Ethr.png");
-    c1->SaveAs("results/far_forward/zdc/zdc_neutron_Ethr.pdf");
-  }
-  {
-    TCanvas* c2 = new TCanvas("c2", "c2");
-    auto h = hThetathr->DrawCopy();
-    h->SetLineWidth(3);
-    h->SetLineColor(kBlue);
-    c2->SaveAs("results/far_forward/zdc/zdc_neutron_Thetathr.png");
-    c2->SaveAs("results/far_forward/zdc/zdc_neutron_Thetathr.pdf");
-  }
-  {
-    TCanvas* c3 = new TCanvas("c3", "c3");
-    auto h = hPhithr->DrawCopy();
-    h->SetLineWidth(3);
-    h->SetLineColor(kBlue);
-    c3->SaveAs("results/far_forward/zdc/zdc_neutron_Phithr.png");
-    c3->SaveAs("results/far_forward/zdc/zdc_neutron_Phithr.pdf");
-  }
-  {
-    TCanvas* c4 = new TCanvas("c4", "c4");
-    auto h = hEtathr->DrawCopy();
-    h->SetLineWidth(3);
-    h->SetLineColor(kBlue);
-    c4->SaveAs("results/far_forward/zdc/zdc_neutron_Etathr.png");
-    c4->SaveAs("results/far_forward/zdc/zdc_neutron_Etathr.pdf");
-  }
-  {
-    TCanvas* c5 = new TCanvas("c5", "c5");
-    auto h = hPhiThetathr->DrawCopy("COLZ");
-    c5->SaveAs("results/far_forward/zdc/zdc_neutron_PhiThetathr.png");
-    c5->SaveAs("results/far_forward/zdc/zdc_neutron_PhiThetathr.pdf");
-  } 
-}
-
diff --git a/benchmarks/far_forward/analysis/analysis_zdc_particles.cxx b/benchmarks/far_forward/analysis/analysis_zdc_particles.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..516f9a00c07726c649f445a84b10f7111248602b
--- /dev/null
+++ b/benchmarks/far_forward/analysis/analysis_zdc_particles.cxx
@@ -0,0 +1,306 @@
+////////////////////////////////////////
+// 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 "eicd/ClusterCollection.h"
+#include "eicd/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(libeicd.so)
+#include "fmt/core.h"
+#include "DD4hep/Detector.h"
+#include "DDG4/Geant4Data.h"
+#include "DDRec/CellIDPositionConverter.h"
+
+using ROOT::RDataFrame;
+
+void analysis_zdc_particles(
+  const std::string& input_fname = "sim_zdc_uniform.edm4hep.root",
+  const std::string& results_path = "results/far_forward/zdc/"
+) {
+  // 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<eicd::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<eicd::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<eicd::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",      {"ZDCHcalClusters.energy"})
+    .Define("Theta_Hcal",  Thetarec, {"ZDCHcalClusters"})
+    .Define("Phi_Hcal",    Phirec,   {"ZDCHcalClusters"})
+	  .Define("Eta_Hcal",    Etarec,   {"ZDCHcalClusters"})
+	;
+
+  // 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/analysis_zdc_photons.cxx b/benchmarks/far_forward/analysis/analysis_zdc_photons.cxx
deleted file mode 100644
index b86527de9de51e6247fc694ad0de431956371429..0000000000000000000000000000000000000000
--- a/benchmarks/far_forward/analysis/analysis_zdc_photons.cxx
+++ /dev/null
@@ -1,156 +0,0 @@
-////////////////////////////////////////
-// Read ROOT output file
-// Plot variables
-// Jihee Kim 09/2021
-////////////////////////////////////////
-
-#include "ROOT/RDataFrame.hxx"
-#include <iostream>
-#include <vector>
-
-#include "edm4hep/MCParticleCollection.h"
-#include "edm4hep/CalorimeterHitCollection.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(libeicd.so)
-#include "fmt/core.h"
-#include "DD4hep/Detector.h"
-#include "DDG4/Geant4Data.h"
-#include "DDRec/CellIDPositionConverter.h"
-
-using ROOT::RDataFrame;
-
-void analysis_zdc_photons(const char* input_fname = "sim_zdc_uniform_photon.edm4hep.root")
-{
-  // 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 Ethr = [](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 [mrad]
-  auto Thetathr = [](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*1000.0;
-  };
-
-  // Phi [rad]
-  auto Phithr = [](const std::vector<edm4hep::MCParticleData> & input) {
-    auto p = input[2];
-    auto phi = std::atan2(p.momentum.y, p.momentum.x);
-    return phi;
-  };
-  
-  // Eta
-  auto Etathr = [](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;
-  };  
-
-  // Define variables
-  auto d1 = d0.Define("Ethr",      Ethr,     {"MCParticles"})
-              .Define("Thetathr",  Thetathr, {"MCParticles"})
-              .Define("Phithr",    Phithr,   {"MCParticles"})
-	      .Define("Etathr",    Etathr,   {"MCParticles"})
-	      ;
-
-  // Define Histograms
-  auto hEthr        = d1.Histo1D({"hEthr",        "Thrown Energy; E_{thr} [GeV]; Events",                100, 0.0, 0.3},  "Ethr");
-  auto hThetathr    = d1.Histo1D({"hThetathr",    "Thrown Theta; #theta_{thr} [mrad]; Events",           100, 0.0, 5.0},  "Thetathr");
-  auto hPhithr      = d1.Histo1D({"hPhithr",      "Thrown Phi; #phi_{thr} [rad]; Events",                100, 0.0, 5.0},  "Phithr");
-  auto hEtathr      = d1.Histo1D({"hEtathr",      "Thrown Eta; #eta_{thr}; Events",                      100, 6.0, 15.0}, "Etathr");
-  auto hPhiThetathr = d1.Histo2D({"hPhiThetathr", "Thrown Phi vs Theta; #theta_{thr} [mrad]; #phi_{thr}", 20, 0.0, 5.0, 20, 0.0, 5.0}, 
-                  "Thetathr", "Phithr");
-
-  // 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 = hEthr->DrawCopy();
-    h->SetLineWidth(3);
-    h->SetLineColor(kBlue);
-    c1->SaveAs("results/far_forward/zdc/zdc_photon_Ethr.png");
-    c1->SaveAs("results/far_forward/zdc/zdc_photon_Ethr.pdf");
-  }
-  {
-    TCanvas* c2 = new TCanvas("c2", "c2");
-    auto h = hThetathr->DrawCopy();
-    h->SetLineWidth(3);
-    h->SetLineColor(kBlue);
-    c2->SaveAs("results/far_forward/zdc/zdc_photon_Thetathr.png");
-    c2->SaveAs("results/far_forward/zdc/zdc_photon_Thetathr.pdf");
-  }
-  {
-    TCanvas* c3 = new TCanvas("c3", "c3");
-    auto h = hPhithr->DrawCopy();
-    h->SetLineWidth(3);
-    h->SetLineColor(kBlue);
-    c3->SaveAs("results/far_forward/zdc/zdc_photon_Phithr.png");
-    c3->SaveAs("results/far_forward/zdc/zdc_photon_Phithr.pdf");
-  }
-  {
-    TCanvas* c4 = new TCanvas("c4", "c4");
-    auto h = hEtathr->DrawCopy();
-    h->SetLineWidth(3);
-    h->SetLineColor(kBlue);
-    c4->SaveAs("results/far_forward/zdc/zdc_photon_Etathr.png");
-    c4->SaveAs("results/far_forward/zdc/zdc_photon_Etathr.pdf");
-  }
-  {
-    TCanvas* c5 = new TCanvas("c5", "c5");
-    auto h = hPhiThetathr->DrawCopy("COLZ");
-    c5->SaveAs("results/far_forward/zdc/zdc_photon_PhiThetathr.png");
-    c5->SaveAs("results/far_forward/zdc/zdc_photon_PhiThetathr.pdf");
-  } 
-}
-
diff --git a/benchmarks/far_forward/analysis/gen_zdc_neutrons.cxx b/benchmarks/far_forward/analysis/gen_zdc_particles.cxx
similarity index 59%
rename from benchmarks/far_forward/analysis/gen_zdc_neutrons.cxx
rename to benchmarks/far_forward/analysis/gen_zdc_particles.cxx
index a755ed4183f37e8a5c9938749a016cc2aa3d6e06..665c38e7b1b6102be5d166d9a2b4a07b3590902a 100644
--- a/benchmarks/far_forward/analysis/gen_zdc_neutrons.cxx
+++ b/benchmarks/far_forward/analysis/gen_zdc_particles.cxx
@@ -1,6 +1,6 @@
 //////////////////////////////////////////////////////////////
 // ZDC detector
-// Single Neutron dataset
+// Single Particle dataset
 // J. Kim 09/2021
 //////////////////////////////////////////////////////////////
 #include "HepMC3/GenEvent.h"
@@ -14,10 +14,12 @@
 #include <math.h>
 #include <random>
 #include <TRandom.h>
+#include <Math/Vector3D.h>
+#include <Math/RotationY.h>
 
 using namespace HepMC3;
 
-void gen_zdc_neutrons(int n_events = 1e6, double e_start = 20.0, double e_end = 140.0, const char* out_fname = "zdc_neutron.hepmc") {
+void gen_zdc_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);
@@ -25,9 +27,12 @@ void gen_zdc_neutrons(int n_events = 1e6, double e_start = 20.0, double e_end =
   // 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.004;
+  double theta_max = 0.015;
 
   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
     // FourVector(px,py,pz,e,pdgid,status)
@@ -38,22 +43,34 @@ void gen_zdc_neutrons(int n_events = 1e6, double e_start = 20.0, double e_end =
     GenParticlePtr p2 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
 
     // Define momentum
-    Double_t p        = r1->Uniform(e_start, e_end);
+    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 px       = p * std::cos(phi) * std::sin(theta);
-    Double_t py       = p * std::sin(phi) * std::sin(theta);
-    Double_t pz       = p * std::cos(theta);
+    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 -25 mrad when particles are being generated
-    //ROOT::Math::XYZVector p0 = {px,py,pz};
-    //ROOT::Math::RotationY r(-0.025);
-    //auto p_rot = r*p0;
+    // 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 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.939565 * 0.939565))), 2112, 1);
+    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 {
+      std::cout << "Particle type " << particle << " not recognized!" << std::endl;
+      exit(-1);
+    }
 
     // Create vertex
     GenVertexPtr v1 = std::make_shared<GenVertex>();
diff --git a/benchmarks/far_forward/analysis/gen_zdc_photons.cxx b/benchmarks/far_forward/analysis/gen_zdc_photons.cxx
deleted file mode 100644
index 753fbae832d3c174cbbabc77a845a4bd407b465e..0000000000000000000000000000000000000000
--- a/benchmarks/far_forward/analysis/gen_zdc_photons.cxx
+++ /dev/null
@@ -1,78 +0,0 @@
-//////////////////////////////////////////////////////////////
-// ZDC detector
-// Single Neutron 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>
-
-using namespace HepMC3;
-
-void gen_zdc_photons(int n_events = 1e6, double e_start = 20.0, double e_end = 140.0, const char* out_fname = "zdc_photon.hepmc") {
-  WriterAscii hepmc_output(out_fname);
-  int events_parsed = 0;
-  GenEvent evt(Units::GEV, Units::MM);
-
-  // Random number generator
-  TRandom* r1 = new TRandom();
-
-  // Set polar angle range [rad]
-  double theta_min = 0.0;
-  double theta_max = 0.004;
-
-  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        = r1->Uniform(e_start, e_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 px       = p * std::cos(phi) * std::sin(theta);
-    Double_t py       = p * std::sin(phi) * std::sin(theta);
-    Double_t pz       = p * std::cos(theta);
-
-    // Rotate the vector in Y by crossing angle -25 mrad when particles are being generated
-    //ROOT::Math::XYZVector p0 = {px,py,pz};
-    //ROOT::Math::RotationY r(-0.025);
-    //auto p_rot = r*p0;
-
-    // FourVector(px,py,pz,e,pdgid,status)
-    // status - type 1 is final state
-    // pdgid 22 - photon
-    GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p)), 22, 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/config.yml b/benchmarks/far_forward/config.yml
index b965eaf677bbe126a2ee17d8eff6a5c320141a02..696ffa693a2fff6fbbb6b36fc41c22939c18b700 100644
--- a/benchmarks/far_forward/config.yml
+++ b/benchmarks/far_forward/config.yml
@@ -11,16 +11,20 @@ B0_far_forward_protons:
   script:
     - bash benchmarks/far_forward/far_forward_protons.sh
 
-ZDC_far_forward_neutrons:
+ZDC_far_forward:
   extends: .rec_benchmark
   stage: run
   needs: ["far_forward:compile"]
   script:
-    - bash benchmarks/far_forward/run_zdc_neutrons.sh
+    - bash benchmarks/far_forward/run_zdc.sh --particle $PARTICLE
+  parallel:
+    matrix:
+      - PARTICLE: ["neutron", "photon"]
 
-ZDC_far_forward_photons:
-  extends: .rec_benchmark
-  stage: run
-  needs: ["far_forward:compile"]
+far_forward:collect:
+  stage: collect
+  needs:
+    - ["B0_far_forward_protons", "ZDC_far_forward"]
   script:
-    - bash benchmarks/far_forward/run_zdc_photons.sh
+    - echo "Done collecting artifacts."
+
diff --git a/benchmarks/far_forward/run_zdc_neutrons.sh b/benchmarks/far_forward/run_zdc.sh
similarity index 82%
rename from benchmarks/far_forward/run_zdc_neutrons.sh
rename to benchmarks/far_forward/run_zdc.sh
index 13b3d6ad1fbe1e2993154154fd953cc2bd0e3ab1..816319e37ce033b991e7d138d9647e64a93f3dce 100644
--- a/benchmarks/far_forward/run_zdc_neutrons.sh
+++ b/benchmarks/far_forward/run_zdc.sh
@@ -1,12 +1,14 @@
 #!/bin/bash
 
 function print_the_help {
-  echo "USAGE: ${0} [--rec-only]  "
+  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=
 
@@ -28,6 +30,11 @@ do
       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"
@@ -69,19 +76,22 @@ fi
 
 export DETECTOR_COMPACT_PATH=${DETECTOR_PATH}/${JUGGLER_DETECTOR_CONFIG}.xml
 
-export JUGGLER_FILE_NAME_TAG="zdc_neutrons"
+export JUGGLER_FILE_NAME_TAG="zdc_${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/${PARTICLE}"
+mkdir -p ${RESULTS_PATH}
+
 echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
 echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
 
 if [[ -z "${REC_ONLY}" && -z "${ANALYSIS_ONLY}" ]] ;
 then
   echo "Generating input events"
-  root -b -q "benchmarks/far_forward/analysis/gen_zdc_neutrons.cxx(${JUGGLER_N_EVENTS}, ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+  root -b -q "benchmarks/far_forward/analysis/gen_zdc_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
@@ -113,11 +123,10 @@ then
   fi
 fi
 
-mkdir -p results/far_forward
-mkdir -p results/far_forward/zdc
+rootls -t ${JUGGLER_REC_FILE}
 
 echo "Running analysis root scripts"
-root -b -q "benchmarks/far_forward/analysis/analysis_zdc_neutrons.cxx+(\"${JUGGLER_REC_FILE}\")"
+root -b -q "benchmarks/far_forward/analysis/analysis_zdc_particles.cxx+(\"${JUGGLER_REC_FILE}\", \"${RESULTS_PATH}\")"
 if [[ "$?" -ne "0" ]] ; then
   echo "ERROR running root analysis script"
   exit 1
@@ -131,4 +140,3 @@ if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
     cp ${JUGGLER_REC_FILE} results/far_forward/zdc/.
   fi
 fi
-
diff --git a/benchmarks/far_forward/run_zdc_photons.sh b/benchmarks/far_forward/run_zdc_photons.sh
deleted file mode 100644
index 9f2d2e20126848ba2eceb3fc6766db77d925839a..0000000000000000000000000000000000000000
--- a/benchmarks/far_forward/run_zdc_photons.sh
+++ /dev/null
@@ -1,134 +0,0 @@
-#!/bin/bash
-
-function print_the_help {
-  echo "USAGE: ${0} [--rec-only]  "
-  echo "  OPTIONS: "
-  echo "            --rec-only   only run gaudi reconstruction part"
-  exit 
-}
-
-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
-      ;;
-    *)    # 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:
-## - JUGGLER_INSTALL_PREFIX:   Install prefix for Juggler (simu/recon)
-## - JUGGLER_DETECTOR:         the detector package we want to use for this benchmark
-## - JUGGLER_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=0.05
-fi
-
-if [[ ! -n  "${E_end}" ]] ; then
-  export E_end=0.25
-fi
-
-if [[ ! -n "${ZDC_SCI_SAMP_FRAC}" ]] ; then
-  export ZDC_PbSCI_SAMP_FRAC=1.0
-fi
-
-export DETECTOR_COMPACT_PATH=${DETECTOR_PATH}/${JUGGLER_DETECTOR_CONFIG}.xml
-
-export JUGGLER_FILE_NAME_TAG="zdc_photons"
-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"
-
-echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
-echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
-
-if [[ -z "${REC_ONLY}" && -z "${ANALYSIS_ONLY}" ]] ;
-then
-  echo "Generating input events"
-  root -b -q "benchmarks/far_forward/analysis/gen_zdc_photons.cxx(${JUGGLER_N_EVENTS}, ${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"
-  ddsim --runType batch \
-    --part.minimalKineticEnergy 0.5*MeV  \
-    --filter.tracker edep0 \
-    -v WARNING \
-    --numberOfEvents ${JUGGLER_N_EVENTS} \
-    --compactFile ${DETECTOR_PATH}/${JUGGLER_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_reconstruction.py
-  if [[ "$?" -ne "0" ]] ; then
-    echo "ERROR running juggler"
-    exit 1
-  fi
-fi
-
-mkdir -p results/far_forward
-mkdir -p results/far_forward/zdc
-
-echo "Running analysis root scripts"
-root -b -q "benchmarks/far_forward/analysis/analysis_zdc_photons.cxx+(\"${JUGGLER_REC_FILE}\")"
-if [[ "$?" -ne "0" ]] ; then
-  echo "ERROR running root analysis script"
-  exit 1
-fi
-
-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/.
-    cp ${JUGGLER_REC_FILE} results/far_forward/zdc/.
-  fi
-fi
-