From 66ac5bf40671931bc80c583ba1463eac6f14767d Mon Sep 17 00:00:00 2001
From: Jihee Kim <jihee.kim@anl.gov>
Date: Mon, 7 Feb 2022 18:58:24 +0000
Subject: [PATCH] Resolve "Add ZDC benchmark"

---
 .../analysis/analysis_zdc_neutrons.cxx        | 156 ++++++++++++++++++
 .../analysis/analysis_zdc_photons.cxx         | 156 ++++++++++++++++++
 .../far_forward/analysis/gen_zdc_neutrons.cxx |  78 +++++++++
 .../far_forward/analysis/gen_zdc_photons.cxx  |  78 +++++++++
 benchmarks/far_forward/config.yml             |  16 ++
 .../far_forward/options/zdc_reconstruction.py | 110 ++++++++++++
 benchmarks/far_forward/run_zdc_neutrons.sh    | 133 +++++++++++++++
 benchmarks/far_forward/run_zdc_photons.sh     | 133 +++++++++++++++
 8 files changed, 860 insertions(+)
 create mode 100644 benchmarks/far_forward/analysis/analysis_zdc_neutrons.cxx
 create mode 100644 benchmarks/far_forward/analysis/analysis_zdc_photons.cxx
 create mode 100644 benchmarks/far_forward/analysis/gen_zdc_neutrons.cxx
 create mode 100644 benchmarks/far_forward/analysis/gen_zdc_photons.cxx
 create mode 100644 benchmarks/far_forward/options/zdc_reconstruction.py
 create mode 100644 benchmarks/far_forward/run_zdc_neutrons.sh
 create mode 100644 benchmarks/far_forward/run_zdc_photons.sh

diff --git a/benchmarks/far_forward/analysis/analysis_zdc_neutrons.cxx b/benchmarks/far_forward/analysis/analysis_zdc_neutrons.cxx
new file mode 100644
index 00000000..125135b0
--- /dev/null
+++ b/benchmarks/far_forward/analysis/analysis_zdc_neutrons.cxx
@@ -0,0 +1,156 @@
+////////////////////////////////////////
+// Read ROOT output file
+// Plot variables
+// Jihee Kim 09/2021
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+#include <vector>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/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)
+R__LOAD_LIBRARY(libDD4pod.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.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<dd4pod::Geant4ParticleData> & input) {
+    auto p = input[2];
+    auto energy = TMath::Sqrt(p.ps.x * p.ps.x + p.ps.y * p.ps.y + p.ps.z * p.ps.z + p.mass * p.mass);
+    return energy;
+  };
+
+  // Theta [mrad]
+  auto Thetathr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+    auto p = input[2];
+    auto theta = p.ps.theta();
+    return theta*1000.0;
+  };
+
+  // Phi [rad]
+  auto Phithr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+    auto p = input[2];
+    auto phi = p.ps.phi();
+    return phi;
+  };
+  
+  // Eta
+  auto Etathr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+    auto p = input[2];
+    auto eta = p.ps.eta();
+    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_photons.cxx b/benchmarks/far_forward/analysis/analysis_zdc_photons.cxx
new file mode 100644
index 00000000..e0cf3102
--- /dev/null
+++ b/benchmarks/far_forward/analysis/analysis_zdc_photons.cxx
@@ -0,0 +1,156 @@
+////////////////////////////////////////
+// Read ROOT output file
+// Plot variables
+// Jihee Kim 09/2021
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+#include <vector>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/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)
+R__LOAD_LIBRARY(libDD4pod.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.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<dd4pod::Geant4ParticleData> & input) {
+    auto p = input[2];
+    auto energy = TMath::Sqrt(p.ps.x * p.ps.x + p.ps.y * p.ps.y + p.ps.z * p.ps.z + p.mass * p.mass);
+    return energy;
+  };
+
+  // Theta [mrad]
+  auto Thetathr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+    auto p = input[2];
+    auto theta = p.ps.theta();
+    return theta*1000.0;
+  };
+
+  // Phi [rad]
+  auto Phithr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+    auto p = input[2];
+    auto phi = p.ps.phi();
+    return phi;
+  };
+  
+  // Eta
+  auto Etathr = [](const std::vector<dd4pod::Geant4ParticleData> & input) {
+    auto p = input[2];
+    auto eta = p.ps.eta();
+    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_neutrons.cxx
new file mode 100644
index 00000000..a755ed41
--- /dev/null
+++ b/benchmarks/far_forward/analysis/gen_zdc_neutrons.cxx
@@ -0,0 +1,78 @@
+//////////////////////////////////////////////////////////////
+// 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_neutrons(int n_events = 1e6, double e_start = 20.0, double e_end = 140.0, const char* 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 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 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);
+
+    // 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/gen_zdc_photons.cxx b/benchmarks/far_forward/analysis/gen_zdc_photons.cxx
new file mode 100644
index 00000000..753fbae8
--- /dev/null
+++ b/benchmarks/far_forward/analysis/gen_zdc_photons.cxx
@@ -0,0 +1,78 @@
+//////////////////////////////////////////////////////////////
+// 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 3aa931fe..81dda471 100644
--- a/benchmarks/far_forward/config.yml
+++ b/benchmarks/far_forward/config.yml
@@ -5,3 +5,19 @@ B0_far_forward_protons:
   script:
     - compile_analyses.py far_forward
     - bash benchmarks/far_forward/far_forward_protons.sh
+
+ZDC_far_forward_neutrons:
+  extends: .rec_benchmark
+  stage: run
+  timeout: 24 hours
+  script:
+    - compile_analyses.py far_forward
+    - bash benchmarks/far_forward/run_zdc_neutrons.sh
+
+ZDC_far_forward_photons:
+  extends: .rec_benchmark
+  stage: run
+  timeout: 24 hours
+  script:
+    - compile_analyses.py far_forward
+    - bash benchmarks/far_forward/run_zdc_photons.sh
diff --git a/benchmarks/far_forward/options/zdc_reconstruction.py b/benchmarks/far_forward/options/zdc_reconstruction.py
new file mode 100644
index 00000000..db98c294
--- /dev/null
+++ b/benchmarks/far_forward/options/zdc_reconstruction.py
@@ -0,0 +1,110 @@
+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__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
+from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
+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', '../atheta/zdc_photons.root')
+kwargs['output'] = os.environ.get('JUGGLER_REC_FILE', 'rec_zdc_photons.root')
+kwargs['compact'] = os.environ.get('DETECTOR_COMPACT_PATH', '../atheta/athena.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',
+    'ZDCEcalHits',
+    'ZDCHcalHits'
+]
+
+podin = PodioInput('PodioReader', collections=sim_colls, OutputLevel=DEBUG)
+podout = PodioOutput('out', filename=kwargs['output'])
+
+# Digitization
+WSciFi_digi = CalHitDigi('WSciFi_digi',
+        inputHitCollection='ZDCEcalHits',
+        outputHitCollection='ZDCEcalRawHits')
+
+PbSci_digi = CalHitDigi('PbSci_digi',
+        inputHitCollection='ZDCHcalHits',
+        outputHitCollection='ZDCHcalRawHits')
+
+# Reconstruction
+WSciFi_reco = CalHitReco('WSciFi_reco',
+        inputHitCollection=WSciFi_digi.outputHitCollection,
+        outputHitCollection='ZDCEcalRecHits',
+        readoutClass='ZDCEcalHits',
+        localDetFields=['system'])
+
+PbSci_reco = CalHitReco('PbSci_reco',
+        inputHitCollection=PbSci_digi.outputHitCollection,
+        outputHitCollection='ZDCHcalRecHits',
+        readoutClass='ZDCHcalHits',
+        localDetFields=['system'])
+
+# Clustering
+WSciFi_cl = IslandCluster('WSciFi_cl',
+        inputHitCollection=WSciFi_reco.outputHitCollection,
+        outputProtoClusterCollection='ZDCEcalProtoClusters',
+        minClusterCenterEdep=3.*MeV,
+        minClusterHitEdep=0.1*MeV,
+        localDistXY=[50*mm, 50*mm],
+        splitCluster=True)
+
+WSciFi_clreco = RecoCoG('WSciFi_clreco',
+        inputHitCollection=WSciFi_cl.inputHitCollection,
+        inputProtoClusterCollection = WSciFi_cl.outputProtoClusterCollection,
+        outputClusterCollection='ZDCEcalClusters',
+        mcHits="ZDCEcalHits",
+        logWeightBase=3.6)
+
+PbSci_cl = IslandCluster('PbSci_cl',
+        inputHitCollection=PbSci_reco.outputHitCollection,
+        outputProtoClusterCollection='ZDCHcalProtoClusters',
+        minClusterCenterEdep=1.*MeV,
+        minClusterHitEdep=0.1*MeV,
+        localDistXY=[200*mm, 200*mm],
+        splitCluster=False)
+
+PbSci_clreco = RecoCoG('PbSci_clreco',
+        inputHitCollection=PbSci_cl.inputHitCollection,
+        inputProtoClusterCollection = PbSci_cl.outputProtoClusterCollection,
+        outputClusterCollection='ZDCHcalClusters',
+        mcHits="ZDCHcalHits",
+        logWeightBase=3.6,
+        samplingFraction=kwargs['PbSci_sf'])
+
+podout.outputCommands = ['keep *']
+
+ApplicationMgr(
+    TopAlg=[podin,
+        WSciFi_digi, PbSci_digi, 
+        WSciFi_reco, PbSci_reco,
+        WSciFi_cl, WSciFi_clreco, PbSci_cl, PbSci_clreco,
+        podout],
+    EvtSel='NONE',
+    EvtMax=kwargs['nev'],
+    ExtSvc=[podioevent],
+    OutputLevel=DEBUG
+)
+
diff --git a/benchmarks/far_forward/run_zdc_neutrons.sh b/benchmarks/far_forward/run_zdc_neutrons.sh
new file mode 100644
index 00000000..c4072140
--- /dev/null
+++ b/benchmarks/far_forward/run_zdc_neutrons.sh
@@ -0,0 +1,133 @@
+#!/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=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_COMPACT_PATH=${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml
+
+export JUGGLER_FILE_NAME_TAG="zdc_neutrons"
+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}"
+
+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\")"
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running script: generating input events"
+    exit 1
+  fi
+
+  echo "Running Geant4 simulation"
+  npsim --runType batch \
+    --part.minimalKineticEnergy 0.5*MeV  \
+    -v WARNING \
+    --numberOfEvents ${JUGGLER_N_EVENTS} \
+    --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.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_neutrons.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
+
diff --git a/benchmarks/far_forward/run_zdc_photons.sh b/benchmarks/far_forward/run_zdc_photons.sh
new file mode 100644
index 00000000..0baf77ec
--- /dev/null
+++ b/benchmarks/far_forward/run_zdc_photons.sh
@@ -0,0 +1,133 @@
+#!/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}.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}.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"
+  npsim --runType batch \
+    --part.minimalKineticEnergy 0.5*MeV  \
+    -v WARNING \
+    --numberOfEvents ${JUGGLER_N_EVENTS} \
+    --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.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
+
-- 
GitLab