diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 7d6a9b2bcf7661f03f9eff6fabb1ec4801823a51..eb6ca7ae0fb1427d8f07a0374e235bd8eaecb45c 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -80,6 +80,7 @@ get_data:
 
 include: 
   - local: 'benchmarks/barrel_ecal/config.yml'
+  - local: 'benchmarks/barrel_hcal/config.yml'
   - local: 'benchmarks/roman_pots/config.yml'
   - local: 'benchmarks/zdc/config.yml'
   - local: 'benchmarks/crystal_calorimeter/config.yml'
@@ -88,7 +89,7 @@ include:
 deploy_results:
   stage: deploy
   needs:
-    - ["collect_results:zdc","collect_results:barrel_ecal","collect_results:crystal_calorimeter"]
+    - ["collect_results:zdc","collect_results:barrel_ecal","collect_results:barrel_hcal","collect_results:crystal_calorimeter"]
   script:
     - echo "deploy results!"
 
diff --git a/benchmarks/barrel_hcal/config.yml b/benchmarks/barrel_hcal/config.yml
new file mode 100644
index 0000000000000000000000000000000000000000..d47f63da9b7f197eff4aefba1598cc7292865375
--- /dev/null
+++ b/benchmarks/barrel_hcal/config.yml
@@ -0,0 +1,107 @@
+sim:hcal_barrel_pions:
+  extends: .det_benchmark
+  stage: simulate
+  script:
+    - bash benchmarks/barrel_hcal/run_hcal_barrel_particles.sh piplus
+    - bash benchmarks/barrel_hcal/run_hcal_barrel_particles.sh piminus
+
+sim:hcal_barrel_kaons:
+  extends: .det_benchmark
+  stage: simulate
+  script:
+    - bash benchmarks/barrel_hcal/run_hcal_barrel_particles.sh kplus
+    - bash benchmarks/barrel_hcal/run_hcal_barrel_particles.sh kminus
+    - bash benchmarks/barrel_hcal/run_hcal_barrel_particles.sh kshort
+    - bash benchmarks/barrel_hcal/run_hcal_barrel_particles.sh klong
+
+sim:hcal_barrel_muons:
+  extends: .det_benchmark
+  stage: simulate
+  script:
+    - bash benchmarks/barrel_hcal/run_hcal_barrel_energy_scan.sh muon
+    - bash benchmarks/barrel_hcal/run_hcal_barrel_particles.sh muon
+
+sim:hcal_barrel_antimuons:
+  extends: .det_benchmark
+  stage: simulate
+  script:
+    - bash benchmarks/barrel_hcal/run_hcal_barrel_energy_scan.sh antimuon
+    - bash benchmarks/barrel_hcal/run_hcal_barrel_particles.sh antimuon
+
+sim:hcal_barrel_protons:
+  extends: .det_benchmark
+  stage: simulate
+  script:
+    - bash benchmarks/barrel_hcal/run_hcal_barrel_energy_scan.sh proton
+    - bash benchmarks/barrel_hcal/run_hcal_barrel_particles.sh proton
+
+bench:hcal_barrel_protons:
+  extends: .det_benchmark
+  stage: benchmarks
+  needs:
+    - ["sim:hcal_barrel_protons"]
+  script:
+    - ls -lhtR sim_output/
+    - rootls -t sim_output/sim_hcal_barrel_proton.root
+    - root -b -q 'benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx+("proton")'
+    - root -b -q 'benchmarks/barrel_hcal/scripts/hcal_barrel_energy_scan_analysis.cxx+("proton")'
+
+bench:hcal_barrel_muons:
+  extends: .det_benchmark
+  stage: benchmarks
+  needs:
+    - ["sim:hcal_barrel_muons"]
+  script:
+    - ls -lhtR sim_output/
+    - rootls -t sim_output/sim_hcal_barrel_muon.root
+    - root -b -q 'benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx+("muon")'
+    - root -b -q 'benchmarks/barrel_hcal/scripts/hcal_barrel_energy_scan_analysis.cxx+("muon")'
+
+bench:hcal_barrel_antimuons:
+  extends: .det_benchmark
+  stage: benchmarks
+  needs:
+    - ["sim:hcal_barrel_antimuons"]
+  script:
+    - ls -lhtR sim_output/
+    - rootls -t sim_output/sim_hcal_barrel_antimuon.root
+    - root -b -q 'benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx+("antimuon")'
+    - root -b -q 'benchmarks/barrel_hcal/scripts/hcal_barrel_energy_scan_analysis.cxx+("antimuon")'
+
+bench:hcal_barrel_kaons:
+  extends: .det_benchmark
+  stage: benchmarks
+  needs:
+    - ["sim:hcal_barrel_kaons"]
+  script:
+    - ls -lhtR sim_output/
+    - rootls -t sim_output/sim_hcal_barrel_kplus.root
+    - root -b -q 'benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx+("kplus")'
+    - rootls -t sim_output/sim_hcal_barrel_kminus.root
+    - root -b -q 'benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx+("kminus")'
+    - rootls -t sim_output/sim_hcal_barrel_kshort.root
+    - root -b -q 'benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx+("kshort")'
+    - rootls -t sim_output/sim_hcal_barrel_klong.root
+    - root -b -q 'benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx+("klong")'
+
+bench:hcal_barrel_pions:
+  extends: .det_benchmark
+  stage: benchmarks
+  needs:
+    - ["sim:hcal_barrel_pions"]
+  script:
+    - ls -lhtR sim_output/
+    - rootls -t sim_output/sim_hcal_barrel_piplus.root
+    - root -b -q 'benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx+("piplus")'
+    - rootls -t sim_output/sim_hcal_barrel_piminus.root
+    - root -b -q 'benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx+("piminus")'
+
+collect_results:barrel_hcal:
+  extends: .det_benchmark
+  stage: collect
+  needs: 
+    - ["bench:hcal_barrel_muons", "bench:hcal_barrel_protons", "bench:hcal_barrel_kaons", "bench:hcal_barrel_pions"]
+  script:
+    - ls -lrht
+    - echo " FIX ME" 
+
diff --git a/benchmarks/barrel_hcal/run_hcal_barrel_energy_scan.sh b/benchmarks/barrel_hcal/run_hcal_barrel_energy_scan.sh
new file mode 100755
index 0000000000000000000000000000000000000000..794333896dbbe35c233b419af962f09102d5a626
--- /dev/null
+++ b/benchmarks/barrel_hcal/run_hcal_barrel_energy_scan.sh
@@ -0,0 +1,24 @@
+#!/bin/bash
+
+export PARTICLE=$1
+E_file="sim_output/hcal_barrel_energy_scan_points_${PARTICLE}.txt"
+
+
+#for E in 1 2 6 10
+for E in 0.25 0.5 1 2 3 4 7 15 20
+do
+   export E_START="$E"
+   export E_END="$E"
+   bash benchmarks/barrel_hcal/run_hcal_barrel_particles.sh "${PARTICLE}" && echo "$E" >> "$E_file" || exit 1
+   path_rootfiles="sim_output/energy_scan/${E}/"
+   path_plots="results/energy_scan/${E}/"
+   mkdir -p "$path_rootfiles"
+   mkdir -p "$path_plots"
+   ls -lthaR sim_output/
+   mv "sim_output/sim_hcal_barrel_${PARTICLE}.root" "$path_rootfiles"
+done
+
+ls -lthaR sim_output
+cat "$E_file"
+
+exit 0
diff --git a/benchmarks/barrel_hcal/run_hcal_barrel_particles.sh b/benchmarks/barrel_hcal/run_hcal_barrel_particles.sh
new file mode 100755
index 0000000000000000000000000000000000000000..edeaba1e4ece7d1cde7dbbf931a3f91d28f43143
--- /dev/null
+++ b/benchmarks/barrel_hcal/run_hcal_barrel_particles.sh
@@ -0,0 +1,66 @@
+#!/bin/bash
+
+if [[ ! -n  "${JUGGLER_DETECTOR}" ]] ; then 
+  export JUGGLER_DETECTOR="topside"
+fi
+
+if [[ ! -n  "${JUGGLER_N_EVENTS}" ]] ; then 
+  export JUGGLER_N_EVENTS=100
+fi
+
+if [[ ! -n  "${E_START}" ]] ; then
+  export E_START=5.0
+fi
+
+if [[ ! -n  "${E_END}" ]] ; then
+  export E_END=5.0
+fi
+
+export PARTICLE=$1
+if [[ ! -n  "${PARTICLE}" ]] ; then
+  export PARTICLE="electron"
+fi
+
+export JUGGLER_FILE_NAME_TAG="hcal_barrel_${PARTICLE}"
+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}"
+
+# Generate the input events
+root -b -q "benchmarks/barrel_hcal/scripts/hcal_barrel_particles_gen.cxx+(${JUGGLER_N_EVENTS}, ${E_START}, ${E_END}, \"${PARTICLE}\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running script: generating input events"
+  exit 1
+fi
+# Plot the input events
+root -b -q "benchmarks/barrel_hcal/scripts/hcal_barrel_particles_reader.cxx+(\"${PARTICLE}\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running script: plotting input events"
+  exit 1
+fi
+
+ls -ltRhL
+
+npsim --runType batch \
+      -v WARNING \
+      --part.minimalKineticEnergy 0.5*GeV  \
+      --numberOfEvents ${JUGGLER_N_EVENTS} \
+      --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
+      --inputFiles data/${JUGGLER_FILE_NAME_TAG}.hepmc \
+      --outputFile sim_output/${JUGGLER_SIM_FILE}
+
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running npdet"
+  exit 1
+fi
+
+# Directory for plots
+mkdir -p results
+
+# Move ROOT output file
+#mv ${JUGGLER_REC_FILE} sim_output/
+
diff --git a/benchmarks/barrel_hcal/scripts/hcal_barrel_common_functions.h b/benchmarks/barrel_hcal/scripts/hcal_barrel_common_functions.h
new file mode 100644
index 0000000000000000000000000000000000000000..3d3efac578ea88865e97b096cb469c468000c0df
--- /dev/null
+++ b/benchmarks/barrel_hcal/scripts/hcal_barrel_common_functions.h
@@ -0,0 +1,23 @@
+//////////////////////////
+// Common Particle Functions
+// M. Scott 05/2021
+//////////////////////////
+
+std::tuple <int, double> extract_particle_parameters(std::string particle_name) {
+    if (particle_name == "electron") return std::make_tuple(11,    0.51099895e-3);
+    if (particle_name == "photon")   return std::make_tuple(22,    0.0);
+    if (particle_name == "positron") return std::make_tuple(-11,   0.51099895e-3);
+    if (particle_name == "proton")   return std::make_tuple(2212,  0.938272);
+    if (particle_name == "muon")     return std::make_tuple(13,    0.1056583745);
+    if (particle_name == "antimuon") return std::make_tuple(-13,   0.1056583745);
+    if (particle_name == "pi0")      return std::make_tuple(111,   0.1349768);
+    if (particle_name == "piplus")   return std::make_tuple(211,   0.13957039);
+    if (particle_name == "piminus")  return std::make_tuple(-211,  0.13957039);
+    if (particle_name == "kplus")    return std::make_tuple(321,   0.493677);
+    if (particle_name == "kminus")   return std::make_tuple(-321,  0.493677);
+    if (particle_name == "kshort")   return std::make_tuple(310,   0.497648);
+    if (particle_name == "klong")    return std::make_tuple(130,   0.497648);
+
+    std::cout << "wrong particle name" << std::endl;
+    abort();
+}
diff --git a/benchmarks/barrel_hcal/scripts/hcal_barrel_energy_scan_analysis.cxx b/benchmarks/barrel_hcal/scripts/hcal_barrel_energy_scan_analysis.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..8ca904ef1193c196c26630111b9ff2f51eee8ff3
--- /dev/null
+++ b/benchmarks/barrel_hcal/scripts/hcal_barrel_energy_scan_analysis.cxx
@@ -0,0 +1,226 @@
+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+#include <fmt/core.h>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/CalorimeterHitCollection.h"
+
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TMath.h"
+#include "TH1.h"
+#include "TF1.h"
+#include "TH1D.h"
+#include "TGraphErrors.h"
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+
+void save_canvas(TCanvas* c, std::string label)
+{
+  c->SaveAs(fmt::format("results/energy_scan/{}.png",label).c_str());
+  c->SaveAs(fmt::format("results/energy_scan/{}.pdf",label).c_str());
+}
+
+void save_canvas(TCanvas* c, std::string label, double E)
+{
+  std::string label_with_E = fmt::format("{}/{}", E, label); 
+  save_canvas(c, label_with_E);
+}
+void save_canvas(TCanvas* c, std::string label, std::string E_label)
+{
+  std::string label_with_E = fmt::format("{}/{}", E_label, label); 
+  save_canvas(c, label_with_E);
+}
+void save_canvas(TCanvas* c, std::string var_label, std::string E_label, std::string particle_label)
+{
+  std::string label_with_E = fmt::format("{}/hcal_barrel_{}_{}", E_label, particle_label, var_label); 
+  save_canvas(c, label_with_E);
+}
+
+std::tuple <double, double, double, double> extract_sampling_fraction_parameters(std::string particle_label, std::string E_label)
+{
+  std::string input_fname = fmt::format("sim_output/energy_scan/{}/sim_hcal_barrel_{}.root", E_label, particle_label);
+  ROOT::EnableImplicitMT();
+  ROOT::RDataFrame d0("events", input_fname);
+
+  // Thrown Energy [GeV]
+  auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+    auto p = input[2];
+    auto energy = TMath::Sqrt(p.psx * p.psx + p.psy * p.psy + p.psz * p.psz + p.mass * p.mass);
+    return energy;
+  };
+
+  // Number of hits
+  auto nhits = [] (const std::vector<dd4pod::CalorimeterHitData>& evt) {return (int) evt.size(); };
+
+  // Energy deposition [GeV]
+  auto Esim = [](const std::vector<dd4pod::CalorimeterHitData>& evt) {
+    auto total_edep = 0.0;
+    for (const auto& i: evt)
+      total_edep += i.energyDeposit;
+    return total_edep;
+  };
+
+  // Sampling fraction = Esampling / Ethrown
+  auto fsam = [](const double sampled, const double thrown) {
+    return sampled / thrown;
+  };
+
+  // Define variables
+  auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
+                .Define("nhits", nhits, {"EcalBarrelHits"})
+                .Define("Esim", Esim, {"EcalBarrelHits"})
+                .Define("fsam", fsam, {"Esim", "Ethr"});
+
+  // Define Histograms
+  auto hEthr = d1.Histo1D(
+      {"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 25.0},
+      "Ethr");
+  auto hNhits =
+      d1.Histo1D({"hNhits", "Number of hits per events; Number of hits; Events",
+                  100, 0.0, 2000.0},
+                 "nhits");
+  auto hEsim = d1.Histo1D(
+      {"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 500, 0.0, 0.5},
+      "Esim");
+  auto hfsam = d1.Histo1D(
+      {"hfsam", "Sampling Fraction; Sampling Fraction; Events", 200, 0.0, 0.05},
+      "fsam");
+
+  // Event Counts
+  auto nevents_thrown = d1.Count();
+  std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
+
+  // Draw Histograms
+  {
+    TCanvas* c1 = new TCanvas("c1", "c1", 700, 500);
+    c1->SetLogy(1);
+    auto h = hEthr->DrawCopy();
+    //h->GetYaxis()->SetTitleOffset(1.4);
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    save_canvas(c1, "Ethr", E_label, particle_label);  
+  }
+
+  {
+    TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
+    c2->SetLogy(1);
+    auto h = hNhits->DrawCopy();
+    //h->GetYaxis()->SetTitleOffset(1.4);
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    save_canvas(c2, "nhits", E_label, particle_label);
+  }
+
+  {
+    TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
+    c3->SetLogy(1);
+    auto h = hEsim->DrawCopy();
+    //h->GetYaxis()->SetTitleOffset(1.4);
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    double up_fit = h->GetMean() + 5*h->GetStdDev();
+    double down_fit = h->GetMean() - 5*h->GetStdDev();
+    h->GetXaxis()->SetRangeUser(0.,up_fit);
+    save_canvas(c3, "Esim", E_label, particle_label);
+  }
+
+  {
+    TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
+    //c4->SetLogy(1);
+    auto h = hfsam->DrawCopy();
+    //h->GetYaxis()->SetTitleOffset(1.4);
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue); 
+    double up_fit = h->GetMean() + 5*h->GetStdDev();
+    double down_fit = h->GetMean() - 5*h->GetStdDev();
+    h->Fit("gaus", "", "", down_fit, up_fit);
+    h->GetXaxis()->SetRangeUser(0.,up_fit);
+    TF1 *gaus = h->GetFunction("gaus");
+    gaus->SetLineWidth(2);
+    gaus->SetLineColor(kRed);    
+    double mean = gaus->GetParameter(1);
+    double sigma = gaus->GetParameter(2);
+    double mean_err = gaus->GetParError(1);
+    double sigma_err = gaus->GetParError(2);
+    save_canvas(c4, "fsam", E_label, particle_label);
+    return std::make_tuple(mean, sigma, mean_err, sigma_err);
+  }
+}
+
+std::vector<std::string> read_scanned_energies(std::string input_energies_fname)
+{
+    std::vector<std::string> scanned_energies;
+    std::string E_label;
+    ifstream E_file (input_energies_fname);
+    if (E_file.is_open())
+    {
+        while (E_file >> E_label)
+        {
+            scanned_energies.push_back(E_label);
+        }
+        E_file.close();
+        return scanned_energies;
+    }
+    else
+    {
+        std::cout << "Unable to open file " << input_energies_fname << std::endl;
+        abort();
+    }
+}
+
+void hcal_barrel_energy_scan_analysis(std::string particle_label = "electron")
+{
+  // Setting for graphs
+  gROOT->SetStyle("Plain");
+  gStyle->SetOptFit(1);
+  gStyle->SetLineWidth(2);
+  gStyle->SetPadTickX(1);
+  gStyle->SetPadTickY(1);
+  gStyle->SetPadGridX(1);
+  gStyle->SetPadGridY(1);
+  gStyle->SetPadLeftMargin(0.14);
+  gStyle->SetPadRightMargin(0.14);
+
+    auto scanned_energies = read_scanned_energies(fmt::format("sim_output/hcal_barrel_energy_scan_points_{}.txt", particle_label));
+
+    TGraphErrors gr_fsam(scanned_energies.size()-1);
+    TGraphErrors gr_fsam_res(scanned_energies.size()-1);
+
+    for (const auto& E_label : scanned_energies) {
+        auto [fsam, fsam_res, fsam_err, fsam_res_err] = extract_sampling_fraction_parameters(particle_label, E_label);
+        auto E = std::stod(E_label);
+
+        gr_fsam.SetPoint(gr_fsam.GetN(),E,100*fsam);
+        gr_fsam.SetPointError(gr_fsam.GetN()-1,0., 100*fsam_err);
+        gr_fsam_res.SetPoint(gr_fsam_res.GetN(),E,100.0*(fsam_res/fsam));
+        auto fsam_res_rel_err = 100.0*(sqrt(pow((fsam_res_err/fsam),2)+pow((fsam_err*fsam_res)/(fsam*fsam),2)));
+        gr_fsam_res.SetPointError(gr_fsam_res.GetN()-1,0.,fsam_res_rel_err);
+    }
+
+    TCanvas* c5 = new TCanvas("c5", "c5", 700, 500);
+    c5->cd();
+    gr_fsam.SetTitle("Sampling Fraction Scan;True Energy [GeV];Sampling Fraction [%]");
+    gr_fsam.SetMarkerStyle(20);
+    gr_fsam.Fit("pol0", "", "", 2., 20.);
+    gr_fsam.Draw("APE");
+    save_canvas(c5, fmt::format("hcal_barrel_{}_fsam_scan", particle_label));
+
+    TCanvas* c6 = new TCanvas("c6", "c6", 700, 500);
+    c6->cd();
+    TF1* func_res = new TF1("func_res", "[0]/sqrt(x) + [1]", 0.25, 20.);
+    func_res->SetLineWidth(2);
+    func_res->SetLineColor(kRed);
+    gr_fsam_res.SetTitle("Energy Resolution;True Energy [GeV];#Delta E/E [%]");
+    gr_fsam_res.SetMarkerStyle(20);
+    gr_fsam_res.Fit(func_res,"R");
+    gr_fsam_res.Draw("APE");
+    save_canvas(c6,fmt::format("hcal_barrel_{}_fsam_scan_res", particle_label));
+}
diff --git a/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx b/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..45559e2ae08de59835cbcc652eb17d847a7f828b
--- /dev/null
+++ b/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx
@@ -0,0 +1,148 @@
+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+#include <fmt/core.h>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/CalorimeterHitCollection.h"
+
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TMath.h"
+#include "TH1.h"
+#include "TF1.h"
+#include "TH1D.h"
+
+#include "hcal_barrel_common_functions.h"
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+void save_canvas(TCanvas* c, std::string label)
+{
+  c->SaveAs(fmt::format("results/{}.png",label).c_str());
+  c->SaveAs(fmt::format("results/{}.pdf",label).c_str());
+}
+
+void save_canvas(TCanvas* c, std::string label, std::string particle_label)
+{
+  std::string label_with_E = fmt::format("hcal_barrel_{}_{}", particle_label, label); 
+  save_canvas(c, label_with_E);
+}
+
+void hcal_barrel_particles_analysis(std::string particle_name = "electron")
+{
+  // Setting for graphs
+  gROOT->SetStyle("Plain");
+  gStyle->SetOptFit(1);
+  gStyle->SetLineWidth(2);
+  gStyle->SetPadTickX(1);
+  gStyle->SetPadTickY(1);
+  gStyle->SetPadGridX(1);
+  gStyle->SetPadGridY(1);
+  gStyle->SetPadLeftMargin(0.14);
+  gStyle->SetPadRightMargin(0.14);
+
+  ROOT::EnableImplicitMT();
+  std::string input_fname = fmt::format("sim_output/sim_hcal_barrel_{}.root", particle_name);
+  ROOT::RDataFrame d0("events", input_fname);
+
+  // Thrown Energy [GeV]
+  auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+    auto p = input[2];
+    auto energy = TMath::Sqrt(p.psx * p.psx + p.psy * p.psy + p.psz * p.psz + p.mass * p.mass);
+    return energy;
+  };
+
+  // Number of hits
+  auto nhits = [] (const std::vector<dd4pod::CalorimeterHitData>& evt) {return (int) evt.size(); };
+
+  // Energy deposition [GeV]
+  auto Esim = [](const std::vector<dd4pod::CalorimeterHitData>& evt) {
+    auto total_edep = 0.0;
+    for (const auto& i: evt)
+      total_edep += i.energyDeposit;
+    return total_edep;
+  };
+
+  // Sampling fraction = Esampling / Ethrown
+  auto fsam = [](const double sampled, const double thrown) {
+    return sampled / thrown;
+  };
+
+  // Define variables
+  auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
+                .Define("nhits", nhits, {"HcalBarrelHits"})
+                .Define("Esim", Esim, {"HcalBarrelHits"})
+                .Define("fsam", fsam, {"Esim", "Ethr"});
+
+  // Define Histograms
+  auto hEthr = d1.Histo1D(
+      {"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 25.0},
+      "Ethr");
+  auto hNhits =
+      d1.Histo1D({"hNhits", "Number of hits per events; Number of hits; Events", 100, 0.0, 2000.0},
+      "nhits");
+  auto hEsim = d1.Histo1D(
+      {"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 500, 0.0, 0.5},
+      "Esim");
+  auto hfsam = d1.Histo1D(
+      {"hfsam", "Sampling Fraction; Sampling Fraction; Events", 200, 0.0, 0.05},
+      "fsam");
+
+  // Event Counts
+  auto nevents_thrown = d1.Count();
+  std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
+
+  // Draw Histograms
+  {
+    TCanvas* c1 = new TCanvas("c1", "c1", 700, 500);
+    c1->SetLogy(1);
+    auto h = hEthr->DrawCopy();
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    save_canvas(c1,"Ethr",particle_name);
+  }
+
+  {
+    TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
+    c2->SetLogy(1);
+    auto h = hNhits->DrawCopy();
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    save_canvas(c2,"nhits",particle_name);
+  }
+
+  {
+    TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
+    c3->SetLogy(1);
+    auto h = hEsim->DrawCopy();
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    double up_fit = h->GetMean() + 5*h->GetStdDev();
+    double down_fit = h->GetMean() - 5*h->GetStdDev();
+    h->GetXaxis()->SetRangeUser(0.,up_fit);
+    save_canvas(c3,"Esim",particle_name);
+  }
+
+  {
+    TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
+    c4->SetLogy(1);
+    auto h = hfsam->DrawCopy();
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    double up_fit = h->GetMean() + 5*h->GetStdDev();
+    double down_fit = h->GetMean() - 5*h->GetStdDev();
+    h->Fit("gaus", "", "", down_fit, up_fit);
+    h->GetXaxis()->SetRangeUser(0.,up_fit);
+    TF1 *gaus = h->GetFunction("gaus");
+    gaus->SetLineWidth(2);
+    gaus->SetLineColor(kRed);
+    save_canvas(c4,"fsam",particle_name);
+  }
+}
+
diff --git a/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_gen.cxx b/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_gen.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..5c9142e08ac5acc48fe3b217db0c403d0145982f
--- /dev/null
+++ b/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_gen.cxx
@@ -0,0 +1,83 @@
+//////////////////////////////////////////////////////////////
+// HCAL Barrel detector
+// Generate particles with particle gun
+// J.Kim 04/02/21
+// M.Zurek 05/05/21
+// W.Deconinck 05/27/21
+//////////////////////////////////////////////////////////////
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/Print.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+
+#include "TMath.h"
+#include "TRandom.h"
+
+#include <cmath>
+#include <iostream>
+#include <math.h>
+#include <random>
+#include <fmt/core.h>
+
+#include "hcal_barrel_common_functions.h"
+
+using namespace HepMC3;
+
+void hcal_barrel_particles_gen(int n_events = 1e6, double e_start = 0.0, double e_end = 20.0, std::string particle_name = "muon") {
+  std::string out_fname = fmt::format("./data/hcal_barrel_{}.hepmc", particle_name);
+  WriterAscii hepmc_output(out_fname);
+  int events_parsed = 0;
+  GenEvent evt(Units::GEV, Units::MM);
+
+  // Random number generator
+  TRandom* r1 = new TRandom();
+
+  // Constraining the solid angle, but larger than that subtended by the
+  // detector
+  // https://indico.bnl.gov/event/7449/contributions/35966/attachments/27177/41430/EIC-DWG-Calo-03192020.pdf
+  // See a figure on slide 26
+  double cos_theta_min = std::cos(M_PI * (45.0 / 180.0));
+  double cos_theta_max = std::cos(M_PI * (135.0 / 180.0));
+
+  for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
+    // FourVector(px,py,pz,e,pdgid,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 costheta = r1->Uniform(cos_theta_min, cos_theta_max);
+    Double_t theta    = std::acos(costheta);
+    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);
+
+    // type 1 is final state
+    auto [id, mass] = extract_particle_parameters(particle_name);
+    GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (mass * mass))), id, 1);
+
+    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 % 10000 == 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/barrel_hcal/scripts/hcal_barrel_particles_reader.cxx b/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_reader.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..6c16a91f29d76c6cd6ecd1998bd92d656f74d4d7
--- /dev/null
+++ b/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_reader.cxx
@@ -0,0 +1,145 @@
+//////////////////////////
+// HCAL Barrel detector
+// J.KIM 04/02/2021
+// M.Zurek 05/05/2021
+// W.Deconinck 05/28/21
+//////////////////////////
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/Print.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+
+#include "TROOT.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TH3.h"
+#include "TStyle.h"
+#include "TCanvas.h"
+#include "TMath.h"
+
+#include <iostream>
+#include <fmt/core.h>
+
+#include "hcal_barrel_common_functions.h"
+
+using namespace HepMC3;
+
+void save_canvas(TCanvas* c, std::string label)
+{
+  c->SaveAs(fmt::format("results/{}.png",label).c_str());
+  c->SaveAs(fmt::format("results/{}.pdf",label).c_str());
+}
+
+void save_canvas(TCanvas* c, std::string label, std::string particle_label)
+{
+  std::string label_with_E = fmt::format("input_hcal_barrel_{}_{}", particle_label, label); 
+  save_canvas(c, label_with_E);
+}
+
+void hcal_barrel_particles_reader(std::string particle_name = "electron")
+{
+  // Setting for graphs
+  gROOT->SetStyle("Plain");
+  gStyle->SetOptFit(1);
+  gStyle->SetLineWidth(1);
+  gStyle->SetPadTickX(1);
+  gStyle->SetPadTickY(1);
+  gStyle->SetPadGridX(1);
+  gStyle->SetPadGridY(1);
+  gStyle->SetPadLeftMargin(0.14);
+  gStyle->SetPadRightMargin(0.17);
+
+  std::string in_fname = fmt::format("./data/hcal_barrel_{}.hepmc",particle_name);
+  ReaderAscii hepmc_input(in_fname);
+  int events_parsed = 0;
+  GenEvent evt(Units::GEV, Units::MM);
+
+  // Histograms
+  TH1F* h_energy = new TH1F(fmt::format("h_{}_energy",particle_name).c_str(), fmt::format("{} energy;E [GeV];Events",particle_name).c_str(),         100, -0.5, 30.5);
+  TH1F* h_eta    = new TH1F(fmt::format("h_{}_eta",particle_name).c_str(),    fmt::format("{} #eta;#eta;Events",particle_name).c_str(),              100, -10.0, 10.0);
+  TH1F* h_theta  = new TH1F(fmt::format("h_{}_theta",particle_name).c_str(),  fmt::format("{} #theta;#theta [degree];Events",particle_name).c_str(), 100, -0.5, 180.5);
+  TH1F* h_phi    = new TH1F(fmt::format("h_{}_phi",particle_name).c_str(),    fmt::format("{} #phi;#phi [degree];Events",particle_name).c_str(),     100, -180.5, 180.5);
+  TH2F* h_pzpt   = new TH2F(fmt::format("h_{}_pzpt",particle_name).c_str(),  fmt::format("{} pt vs pz;pt [GeV];pz [GeV]",particle_name).c_str(),    100, -0.5, 30.5, 100, -30.5, 30.5);
+  TH2F* h_pxpy   = new TH2F(fmt::format("h_{}_pxpy",particle_name).c_str(),  fmt::format("{} px vs py;px [GeV];py [GeV]",particle_name).c_str(),    100, -30.5, 30.5, 100, -30.5, 30.5);
+  TH3F* h_p      = new TH3F(fmt::format("h_{}_p",particle_name).c_str(),     fmt::format("{} p;px [GeV];py [GeV];pz [GeV]",particle_name).c_str(),  100, -30.5, 30.5, 100, -30.5, 30.5, 100, -30.5, 30.5);
+
+  while (!hepmc_input.failed()) {
+    // Read event from input file
+    hepmc_input.read_event(evt);
+    // If reading failed - exit loop
+    if (hepmc_input.failed())
+      break;
+
+    auto [id, mass] = extract_particle_parameters(particle_name);
+
+    for (const auto& v : evt.vertices()) {
+      for (const auto& p : v->particles_out()) {
+        if (p->pid() == id) {
+          h_energy->Fill(p->momentum().e());
+          h_eta->Fill(p->momentum().eta());
+          h_theta->Fill(p->momentum().theta() * TMath::RadToDeg());
+          h_phi->Fill(p->momentum().phi() * TMath::RadToDeg());
+          h_pzpt->Fill(TMath::Sqrt(p->momentum().px() * p->momentum().px() + p->momentum().py() * p->momentum().py()), p->momentum().pz());
+          h_pxpy->Fill(p->momentum().px(), p->momentum().py());
+          h_p->Fill(p->momentum().px(), p->momentum().py(), p->momentum().pz());
+        }
+      }
+    }
+    evt.clear();
+    events_parsed++;
+  }
+  std::cout << "Events parsed and written: " << events_parsed << std::endl;
+
+  TCanvas* c = new TCanvas("c", "c", 500, 500);
+  h_energy->GetYaxis()->SetTitleOffset(1.8);
+  h_energy->SetLineWidth(2);
+  h_energy->SetLineColor(kBlue);
+  h_energy->DrawClone();
+  save_canvas(c, "energy", particle_name);
+
+  TCanvas* c1 = new TCanvas("c1", "c1", 500, 500);
+  h_eta->GetYaxis()->SetTitleOffset(1.9);
+  h_eta->SetLineWidth(2);
+  h_eta->SetLineColor(kBlue);
+  h_eta->DrawClone();
+  save_canvas(c1, "eta", particle_name);
+
+  TCanvas* c2 = new TCanvas("c2", "c2", 500, 500);
+  h_theta->GetYaxis()->SetTitleOffset(1.8);
+  h_theta->SetLineWidth(2);
+  h_theta->SetLineColor(kBlue);
+  h_theta->DrawClone();
+  save_canvas(c2, "theta", particle_name);
+
+  TCanvas* c3 = new TCanvas("c3", "c3", 500, 500);
+  h_phi->GetYaxis()->SetTitleOffset(1.8);
+  h_phi->SetLineWidth(2);
+  h_phi->GetYaxis()->SetRangeUser(0.0, h_phi->GetMaximum() + 100.0);
+  h_phi->SetLineColor(kBlue);
+  h_phi->DrawClone();
+  save_canvas(c3, "phi", particle_name);
+
+  TCanvas* c4 = new TCanvas("c4", "c4", 500, 500);
+  h_pzpt->GetYaxis()->SetTitleOffset(1.4);
+  h_pzpt->SetLineWidth(2);
+  h_pzpt->SetLineColor(kBlue);
+  h_pzpt->DrawClone("COLZ");
+  save_canvas(c4, "pzpt", particle_name);
+
+  TCanvas* c5 = new TCanvas("c5", "c5", 500, 500);
+  h_pxpy->GetYaxis()->SetTitleOffset(1.4);
+  h_pxpy->SetLineWidth(2);
+  h_pxpy->SetLineColor(kBlue);
+  h_pxpy->DrawClone("COLZ");
+  save_canvas(c5, "pxpy", particle_name);
+
+  TCanvas* c6 = new TCanvas("c6", "c6", 500, 500);
+  h_p->GetYaxis()->SetTitleOffset(1.8);
+  h_p->GetXaxis()->SetTitleOffset(1.6);
+  h_p->GetZaxis()->SetTitleOffset(1.6);
+  h_p->SetLineWidth(2);
+  h_p->SetLineColor(kBlue);
+  h_p->DrawClone();
+  save_canvas(c6, "p", particle_name);
+}
+