From c7c10fc04ed3b7419adc7c9269be600963385d44 Mon Sep 17 00:00:00 2001
From: Sylvester Joosten <sjoosten@anl.gov>
Date: Thu, 23 Sep 2021 04:25:03 +0000
Subject: [PATCH] Update for new EICD/Juggler

---
 .../clustering/options/full_cal_reco.py       |  14 +-
 .../clustering/scripts/cluster_plots.py       |   3 -
 benchmarks/ecal/options/barrel.py             |   2 +-
 benchmarks/ecal/options/endcap_e.py           |   2 +-
 benchmarks/ecal/options/endcap_i.py           |   2 +-
 .../ecal/scripts/deprecated/cal_eng_res.C     |  84 ---
 .../deprecated/crystal_cal_electrons.cxx      | 130 -----
 .../deprecated/emcal_barrel_electrons.cxx     |  80 ---
 .../emcal_barrel_electrons_analysis.cxx       | 174 -------
 .../emcal_barrel_electrons_reader.cxx         | 125 -----
 .../scripts/deprecated/emcal_barrel_pions.cxx |  80 ---
 .../emcal_barrel_pions_analysis.cxx           | 174 -------
 .../deprecated/emcal_barrel_pions_reader.cxx  | 125 -----
 .../scripts/deprecated/emcal_electrons.cxx    |  90 ----
 .../deprecated/emcal_electrons_analysis.cxx   | 492 ------------------
 .../deprecated/emcal_electrons_reader.cxx     | 126 -----
 .../ecal/scripts/deprecated/emcal_pi0.cxx     |  87 ----
 .../scripts/deprecated/emcal_pi0_reader.cxx   | 138 -----
 .../deprecated/emcal_pion_anlaysis.cxx        | 251 ---------
 .../deprecated/full_emcal_electrons.cxx       |  89 ----
 .../full_emcal_electrons_analysis.cxx         |  70 ---
 benchmarks/ecal/scripts/deprecated/makeplot.C | 249 ---------
 .../ecal/scripts/deprecated/makeplot_input.C  | 102 ----
 .../ecal/scripts/deprecated/makeplot_pi0.C    | 315 -----------
 .../ecal/scripts/deprecated/rdf_test.cxx      |  66 ---
 benchmarks/ecal/scripts/deprecated/read_eng.C |  70 ---
 .../deprecated/rec_emcal_electrons_reader.C   | 370 -------------
 .../rec_emcal_resolution_analysis.cxx         | 248 ---------
 benchmarks/ecal/scripts/draw_clusters.py      |   3 -
 .../imaging_ecal/options/hybrid_cluster.py    |   4 +-
 .../imaging_ecal/options/imaging_2dcluster.py |   2 +-
 .../options/imaging_topocluster.py            |   2 +-
 .../imaging_ecal/options/scfi_cluster.py      |   2 +-
 .../scripts/aggregator_occupancy.py           |   1 -
 34 files changed, 15 insertions(+), 3757 deletions(-)
 delete mode 100644 benchmarks/ecal/scripts/deprecated/cal_eng_res.C
 delete mode 100644 benchmarks/ecal/scripts/deprecated/crystal_cal_electrons.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons_analysis.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons_reader.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/emcal_barrel_pions.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/emcal_barrel_pions_analysis.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/emcal_barrel_pions_reader.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/emcal_electrons.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/emcal_electrons_analysis.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/emcal_electrons_reader.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/emcal_pi0.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/emcal_pi0_reader.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/emcal_pion_anlaysis.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/full_emcal_electrons.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/full_emcal_electrons_analysis.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/makeplot.C
 delete mode 100644 benchmarks/ecal/scripts/deprecated/makeplot_input.C
 delete mode 100644 benchmarks/ecal/scripts/deprecated/makeplot_pi0.C
 delete mode 100644 benchmarks/ecal/scripts/deprecated/rdf_test.cxx
 delete mode 100644 benchmarks/ecal/scripts/deprecated/read_eng.C
 delete mode 100644 benchmarks/ecal/scripts/deprecated/rec_emcal_electrons_reader.C
 delete mode 100644 benchmarks/ecal/scripts/deprecated/rec_emcal_resolution_analysis.cxx

diff --git a/benchmarks/clustering/options/full_cal_reco.py b/benchmarks/clustering/options/full_cal_reco.py
index a42d4199..cd983a4b 100644
--- a/benchmarks/clustering/options/full_cal_reco.py
+++ b/benchmarks/clustering/options/full_cal_reco.py
@@ -106,7 +106,7 @@ ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
         inputHitCollection=ce_ecal_cl.inputHitCollection,
         inputProtoClusterCollection=ce_ecal_cl.outputProtoClusterCollection,
         outputClusterCollection="EcalEndcapNClusters",
-        outputInfoCollection="EcalEndcapNClustersInfo",
+        mcHits="EcalEndcapNHits",
         samplingFraction=0.998,      # this accounts for a small fraction of leakage
         logWeightBase=4.6)
 
@@ -151,7 +151,7 @@ ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
         inputHitCollection=ci_ecal_cl.inputHitCollection,
         inputProtoClusterCollection=ci_ecal_cl.outputProtoClusterCollection,
         outputClusterCollection="EcalEndcapPClusters",
-        outputInfoCollection="EcalEndcapPClustersInfo",
+        mcHits="EcalEndcapPHits",
         logWeightBase=6.2,
         samplingFraction=ci_ecal_sf)
 
@@ -189,8 +189,8 @@ cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
         samplingFraction=cb_ecal_sf,
         inputHitCollection=cb_ecal_cl.inputHitCollection,
         inputProtoClusterCollection=cb_ecal_cl.outputProtoClusterCollection,
+        mcHits="EcalBarrelHits",
         outputClusterCollection="EcalBarrelImagingClusters",
-        outputInfoCollection="EcalBarrelImagingClustersInfo",
         outputLayerCollection="EcalBarrelImagingLayers")
 
 #Central ECAL SciFi
@@ -237,7 +237,7 @@ scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
        inputHitCollection=scfi_barrel_cl.inputHitCollection,
        inputProtoClusterCollection=scfi_barrel_cl.outputProtoClusterCollection,
        outputClusterCollection="EcalBarrelScFiClusters",
-       outputInfoCollection="EcalBarrelScFiClustersInfo",
+       mcHits="EcalBarrelScFiHits",
        logWeightBase=6.2,
        samplingFraction= scifi_barrel_sf)
 
@@ -281,7 +281,7 @@ cb_hcal_clreco = RecoCoG("cb_hcal_clreco",
         inputHitCollection=cb_hcal_cl.inputHitCollection,
         inputProtoClusterCollection=cb_hcal_cl.outputProtoClusterCollection,
         outputClusterCollection="HcalBarrelClusters",
-        outputInfoCollection="HcalBarrelClustersInfo",
+        mcHits="HcalBarrelHits",
         logWeightBase=6.2,
         samplingFraction=cb_hcal_sf)
 
@@ -321,7 +321,7 @@ ci_hcal_clreco = RecoCoG("ci_hcal_clreco",
         inputHitCollection=ci_hcal_cl.inputHitCollection,
         inputProtoClusterCollection=ci_hcal_cl.outputProtoClusterCollection,
         outputClusterCollection="HcalEndcapPClusters",
-        outputInfoCollection="HcalEndcapPClustersInfo",
+        mcHits="HcalEndcapPHits",
         logWeightBase=6.2,
         samplingFraction=ci_hcal_sf)
 
@@ -361,7 +361,7 @@ ce_hcal_clreco = RecoCoG("ce_hcal_clreco",
         inputHitCollection=ce_hcal_cl.inputHitCollection,
         inputProtoClusterCollection=ce_hcal_cl.outputProtoClusterCollection,
         outputClusterCollection="HcalEndcapNClusters",
-        outputInfoCollection="HcalEndcapNClustersInfo",
+        mcHits="HcalEndcapNHits",
         logWeightBase=6.2,
         samplingFraction=ce_hcal_sf)
 
diff --git a/benchmarks/clustering/scripts/cluster_plots.py b/benchmarks/clustering/scripts/cluster_plots.py
index 29c2ebdd..567c5151 100644
--- a/benchmarks/clustering/scripts/cluster_plots.py
+++ b/benchmarks/clustering/scripts/cluster_plots.py
@@ -25,7 +25,6 @@ def load_root_macros(arg_macros):
 def flatten_collection(rdf, collection, cols=None):
     if not cols:
         cols = [str(c) for c in rdf.GetColumnNames() if str(c).startswith('{}.'.format(collection))]
-        cols += [str(c) for c in rdf.GetColumnNames() if str(c).startswith('{}Info.'.format(collection))]
     else:
         cols = ['{}.{}'.format(collection, c) for c in cols]
     if not cols:
@@ -53,7 +52,6 @@ def thrown_particles_figure(rdf, save, mcbranch="mcparticles"):
     # define truth particle info
     dft = flatten_collection(rdf, mcbranch, ['genStatus', 'pdgID', 'ps.x', 'ps.y', 'ps.z', 'mass'])
     dft.rename(columns={c: c.replace(mcbranch + '.', '') for c in dft.columns}, inplace=True)
-    dft.rename(columns={c: c.replace(mcbranch + 'Info.', '') for c in dft.columns}, inplace=True)
     # select thrown particles
     dft = dft[dft['genStatus'] == 1]
 
@@ -160,7 +158,6 @@ if __name__ == '__main__':
             print('{} do not have valid entries, skip it'.format(coll))
             continue
         df.rename(columns={c: c.replace(coll + '.', '') for c in df.columns}, inplace=True)
-        df.rename(columns={c: c.replace(coll + 'Info.', '') for c in df.columns}, inplace=True)
         # print(df[['eta', 'polar.theta', 'position.x', 'position.y', 'position.z']])
         fig, axs = plt.subplots(2, 2, figsize=(12, 8), dpi=160)
         ncl = df.groupby('event')['ID.value'].nunique().values
diff --git a/benchmarks/ecal/options/barrel.py b/benchmarks/ecal/options/barrel.py
index b78eb191..58b7bef2 100644
--- a/benchmarks/ecal/options/barrel.py
+++ b/benchmarks/ecal/options/barrel.py
@@ -86,7 +86,7 @@ cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
         inputProtoClusterCollection=cb_ecal_cl.outputProtoClusterCollection,
         outputClusterCollection="EcalBarrelImagingClusters",
         outputLayerCollection="EcalBarrelImagingLayers",
-        outputInfoCollection="EcalBarrelImagingClustersInfo")
+        mcHits="EcalBarrelHits")
 
 podout.outputCommands = ['drop *',
         'keep mcparticles2',
diff --git a/benchmarks/ecal/options/endcap_e.py b/benchmarks/ecal/options/endcap_e.py
index 36b9376f..5a29973b 100644
--- a/benchmarks/ecal/options/endcap_e.py
+++ b/benchmarks/ecal/options/endcap_e.py
@@ -78,7 +78,7 @@ ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
         inputHitCollection=ce_ecal_cl.inputHitCollection,
         inputProtoClusterCollection=ce_ecal_cl.outputProtoClusterCollection,
         outputClusterCollection="EcalEndcapNClusters",
-        outputInfoCollection="EcalEndcapNClustersInfo",
+        mcHits="EcalEndcapNHits",
         samplingFraction=0.998,      # this accounts for a small fraction of leakage
         logWeightBase=4.6)
 
diff --git a/benchmarks/ecal/options/endcap_i.py b/benchmarks/ecal/options/endcap_i.py
index 0bd79e2c..1a7fd342 100644
--- a/benchmarks/ecal/options/endcap_i.py
+++ b/benchmarks/ecal/options/endcap_i.py
@@ -87,7 +87,7 @@ ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
         inputHitCollection=ci_ecal_cl.inputHitCollection,
         inputProtoClusterCollection=ci_ecal_cl.outputProtoClusterCollection,
         outputClusterCollection="EcalEndcapPClusters",
-        outputInfoCollection="EcalEndcapPClustersInfo",
+        mcHits="EcalEndcapPHits",
         logWeightBase=6.2,
         samplingFraction=ci_ecal_sf)
 
diff --git a/benchmarks/ecal/scripts/deprecated/cal_eng_res.C b/benchmarks/ecal/scripts/deprecated/cal_eng_res.C
deleted file mode 100644
index 3aa4220a..00000000
--- a/benchmarks/ecal/scripts/deprecated/cal_eng_res.C
+++ /dev/null
@@ -1,84 +0,0 @@
-int cal_eng_res(const char* input_fname = "sim_output/read_eng_output.root")
-{
-  // Setting Figures
-  gROOT->SetStyle("Plain");
-  gStyle->SetLineWidth(3);
-  gStyle->SetOptStat("nem");
-  gStyle->SetOptFit(1);
-  gStyle->SetPadTickX(1);
-  gStyle->SetPadTickY(1);
-  gStyle->SetPadGridX(1);
-  gStyle->SetPadGridY(1);
-  gStyle->SetPadLeftMargin(0.14);
-
-  // Read ROOT File
-  TFile *f = new TFile(input_fname,"r");
-  TTree *tEngRes = (TTree *)f->Get("tEngRes");
-  // Variables from ROOT tree
-  Int_t    ievent;
-  Double_t tot_clust_eng; // total cluster energy [GeV]
-  Int_t    ncluster;      // number of clusters per event
-  Double_t mc_eng;        // thrown energy [GeV]
-  Double_t sim_tru_eng;   // total truth simulated energy [GeV]
-  Double_t sim_eng;       // total simulated energy [GeV]
-  // Read Branches
-  tEngRes->SetBranchAddress("ievent",        &ievent);
-  tEngRes->SetBranchAddress("tot_clust_eng", &tot_clust_eng);
-  tEngRes->SetBranchAddress("ncluster",      &ncluster);
-  tEngRes->SetBranchAddress("mc_eng",        &mc_eng);
-  tEngRes->SetBranchAddress("sim_tru_eng",   &sim_tru_eng);
-  tEngRes->SetBranchAddress("sim_eng",       &sim_eng);
-
-  // Setting for Canvas
-  TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
-  TCanvas *c2 = new TCanvas("c2", "c2", 700, 500);
-
-  // Declare Histrograms
-  TH1D *h1 = new TH1D("hEnergyRes","Energy Resolution", 100,-0.3,0.3);
-  TH1D *h2 = new TH1D("hEnergyResCUT","Energy Resolution with CUT", 100,-0.3,0.3);
-
-  // Variables that used in calculations
-  Double_t diff_energy; 
-  Double_t eng_res;
-
-  Int_t nentries = tEngRes->GetEntries();
-
-  for(Int_t ievent = 0; ievent < nentries; ievent++)
-  {
-	tEngRes->GetEntry(ievent);
-
-	if(ncluster == 1)
-	{
-		diff_energy = tot_clust_eng - mc_eng;
-		eng_res     = (diff_energy / mc_eng);
-		cout << tot_clust_eng << "\t" << mc_eng << "\t" << eng_res << endl;
-		h1->Fill(eng_res, 1.0);
-
-		if(tot_clust_eng > 0.5)
-			h2->Fill(eng_res, 1.0);
-	}
-  }
-
-  // Drawing and Saving Figures
-  c1->cd();
-  h1->SetLineColor(kBlue);
-  h1->SetLineWidth(2);
-  h1->GetXaxis()->SetTitle("#DeltaE/E");
-  h1->GetYaxis()->SetTitle("events");
-  h1->GetYaxis()->SetTitleOffset(1.4);
-  h1->Fit("gaus");
-  h1->Draw();
-  c1->SaveAs("results/hEngRes.png");
-
-  c2->cd();
-  h2->SetLineColor(kBlue);
-  h2->SetLineWidth(2);
-  h2->GetXaxis()->SetTitle("#DeltaE/E");
-  h2->GetYaxis()->SetTitle("events");
-  h2->GetYaxis()->SetTitleOffset(1.4);
-  h2->Fit("gaus");
-  h2->Draw();
-  c2->SaveAs("results/hEngRes_CUT.png");
-
-  return 0;
-}
diff --git a/benchmarks/ecal/scripts/deprecated/crystal_cal_electrons.cxx b/benchmarks/ecal/scripts/deprecated/crystal_cal_electrons.cxx
deleted file mode 100644
index f5d781ea..00000000
--- a/benchmarks/ecal/scripts/deprecated/crystal_cal_electrons.cxx
+++ /dev/null
@@ -1,130 +0,0 @@
-#include "ROOT/RDataFrame.hxx"
-#include <iostream>
-
-#include "dd4pod/Geant4ParticleCollection.h"
-#include "eicd/ClusterCollection.h"
-#include "eicd/ClusterData.h"
-// If root cannot find a dd4pod header make sure to add the include dir to ROOT_INCLUDE_PATH:
-//  export ROOT_INCLUDE_PATH=$HOME/include:$ROOT_INCLUDE_PATH
-
-using ROOT::RDataFrame;
-using namespace ROOT::VecOps;
-
-auto momentum = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
-  std::vector<double> result;
-  for (size_t i = 0; i < in.size(); ++i) {
-   result.push_back(in[i].E());
-  }
-  return result;
-};
-auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
-  std::vector<double> result;
-  for (size_t i = 0; i < in.size(); ++i) {
-   result.push_back(in[i].Theta()*180/M_PI);
-  }
-  return result;
-};
-auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
-  std::vector<ROOT::Math::PxPyPzMVector> result;
-  ROOT::Math::PxPyPzMVector lv;
-  for (size_t i = 0; i < in.size(); ++i) {
-    lv.SetCoordinates(in[i].psx, in[i].psy, in[i].psz, in[i].mass);
-    result.push_back(lv);
-  }
-  return result;
-};
-auto pt = [](std::vector<dd4pod::Geant4ParticleData> const& in) {
-  std::vector<float> result;
-  for (size_t i = 0; i < in.size(); ++i) {
-    result.push_back(std::sqrt(in[i].psx * in[i].psx + in[i].psy * in[i].psy));
-  }
-  return result;
-};
-auto eta = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
-  std::vector<float> result;
-  ROOT::Math::PxPyPzMVector lv;
-  for (size_t i = 0; i < in.size(); ++i) {
-    lv.SetCoordinates(in[i].psx, in[i].psy, in[i].psz, in[i].mass);
-    result.push_back(lv.Eta());
-  }
-  return result;
-};
-
-void crystal_cal_electrons(const char* in_fname = "topside/rec_emcal_uniform_electrons.root")
-{
-  ROOT::EnableImplicitMT();
-  ROOT::RDataFrame df("events", in_fname);
-
-  //TH1D* h1 = new TH1D("h1", "Scattering Angle(#theta)", 100, 130.0, 180.0);
-  //TH1D* h2 = new TH1D("h2", "Pseudo-rapidity(#eta)", 100, -5.0, 0.0);
-  //TH2D* h3 = new TH2D("h3", "Cluster E vs Pseudo-rapidity", 100, e_start - 0.5, e_end + 0.5, 100, -5.0, 0.0);
-  //TH1D* h4 = new TH1D("h4", "Reconstructed energy per event", 100, e_start - 0.5, e_end + 0.5);
-  //TH1D* h5 = new TH1D("h5", "Number of Clusters per event", 5, -0.5, 4.5);
-  //TH1D* h6 = new TH1D("h6", "Scattering Angle(#theta) with CUT", 100, 130.0, 180.0);
-  //TH1D* h7 = new TH1D("h7", "Pseudo-rapidity(#eta) with CUT", 100, -5.0, 0.0);
-  //TH2D *h8  = new TH2D("h8","Cluster Hit Position",                        62,-62.0,62.0,62,-62.0,62.0);
-  //TH2D *h9  = new TH2D("h9","All Hit Position",                            62,-62.0,62.0,62,-62.0,62.0);
-  //TH1D *h10 = new TH1D("hEnergyRes","Energy Resolution",                   100,-0.3,0.3);
-  //TH1D *h11 = new TH1D("hEnergyResCUT","Energy Resolution with CUT",       100,-0.3,0.3);
-  //TH1D *h12 = new TH1D("h12","Thrown momentum",                            61,e_start-0.5,e_end+0.5);
-  //TH1D *h13 = new TH1D("h13","Thrown momentum for reconstructed particle", 61,e_start-0.5,e_end+0.5);
-  //TH1D *h14 = new TH1D("h14","Ratio p_{rec}/p_{thr}",                      61,e_start-0.5,e_end+0.5);
-
-  auto d0 = df.Define("isThrown", "mcparticles2.genStatus == 1")
-                .Define("thrownParticles", "mcparticles2[isThrown]")
-                .Define("thrownP", fourvec, {"thrownParticles"})
-                .Define("thrownEta", eta, {"thrownParticles"})
-                .Define("thrownTheta", theta, {"thrownP"})
-                .Define("thrownMomentum", momentum, {"thrownP"})
-                .Define("Ecluster", [](const std::vector<eic::ClusterData>& in) {
-                  std::vector<float> res;
-                  for (const auto& i : in)
-                    res.push_back(i.energy);
-                  return res;
-                },{"EcalClusters"})
-                .Define("nclusters","EcalClusters.size()") ;
-
-  auto d1 = d0.Filter("nclusters==1");
-  auto c_nclusters1 = d1.Count();
-  auto c_thrown = d0.Count();
-
-  auto h_eta_thrown   = d0.Histo1D({"h_eta_thrown", " ; #eta ", 100, -5.0, 0.0}, "thrownEta");
-  auto h_theta_thrown = d0.Histo1D({"h_theta_thrown", "; #theta", 100, 30.0, 180.0},"thrownTheta");
-  auto h_momentum_thrown = d0.Histo1D({"h_momentum_thrown", "; E [GeV]", 100, 0.0, 30.0},"thrownMomentum");
-  auto h_nclusters = d0.Histo1D({"h_nclusters", "; N clusters", 6, 0,6},"nclusters");
-  auto h_Ecluster = d0.Histo1D({"h_Ecluster", ";  cluster E [GeV]",100, 0,30},"Ecluster");
-
-  auto h_Ecluster1 = d1.Histo1D({"h_Ecluster1", "One cluster events;  cluster E [GeV]",100, 0,30},"Ecluster");
-  auto h_momentum_thrown1 = d1.Histo1D({"h_momentum_thrown", "; E [GeV]", 100, 0.0, 30.0},"thrownMomentum");
-
-  //auto h_Ecluster2 = d1.Filter("thrownMomentum > 4").Histo1D({"h_Ecluster1", "One cluster events;  cluster E [GeV]",100, 0,30},"Ecluster");
-  //auto h_momentum_thrown2 = d1.Filter("thrownMomentum > 4").Histo1D({"h_momentum_thrown", "; E [GeV]", 100, 0.0, 30.0},"thrownMomentum");
-
-  auto c = new TCanvas();
-
-  h_eta_thrown->DrawCopy();
-  c->SaveAs("results/crystal_cal_electrons_etaThrown.png");
-
-  h_theta_thrown->DrawCopy();
-  c->SaveAs("results/crystal_cal_electrons_thetaThrown.png");
-
-  h_nclusters->DrawCopy();
-  c->SaveAs("results/crystal_cal_electrons_nclusters.png");
-
-  h_Ecluster->DrawCopy();
-  h_Ecluster1->SetLineColor(2);
-  h_Ecluster1->DrawCopy("same");
-  //h_Ecluster2->SetLineColor(4);
-  //h_Ecluster2->DrawCopy("same");
-  c->SaveAs("results/crystal_cal_electrons_Ecluster.png");
-
-
-  std::cout <<  (*c_nclusters1) << " / " << (*c_thrown) << " = single cluster events/all \n";
-
-  //c->SaveAs("results/crystal_cal_electrons_Ecluster.png");
-
-
-  //std::string outfilename = "rdf_test.root";
-  //df2.Snapshot("events", outfilename, {"MCParticles_pt", "mcparticles"});
-  return 0;
-}
diff --git a/benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons.cxx b/benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons.cxx
deleted file mode 100644
index 12a78da3..00000000
--- a/benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons.cxx
+++ /dev/null
@@ -1,80 +0,0 @@
-//////////////////////////////////////////////////////////////
-// EMCAL Barrel detector
-// Single Electron dataset
-// J.KIM 04/02/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>
-
-using namespace HepMC3;
-
-void emcal_barrel_electrons(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/emcal_barrel_electrons.hepmc") {
-  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 111 - pi0
-    // 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);
-
-    // Generates random vectors, uniformly distributed over the surface of a
-    // sphere of given radius, in this case momentum.
-    // r1->Sphere(px, py, pz, p);
-
-    // type 1 is final state
-    // pdgid 11 - electron 0.510 MeV/c^2
-    GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.000511 * 0.000511))), 11, 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/ecal/scripts/deprecated/emcal_barrel_electrons_analysis.cxx b/benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons_analysis.cxx
deleted file mode 100644
index bab641eb..00000000
--- a/benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons_analysis.cxx
+++ /dev/null
@@ -1,174 +0,0 @@
-////////////////////////////////////////
-// Read reconstruction ROOT output file
-// Plot variables
-////////////////////////////////////////
-
-#include "ROOT/RDataFrame.hxx"
-#include <iostream>
-
-#include "dd4pod/CalorimeterHitCollection.h"
-#include "dd4pod/Geant4ParticleCollection.h"
-#include "eicd/CalorimeterHitCollection.h"
-#include "eicd/CalorimeterHitData.h"
-#include "eicd/ClusterCollection.h"
-#include "eicd/ClusterData.h"
-#include "eicd/RawCalorimeterHitCollection.h"
-#include "eicd/RawCalorimeterHitData.h"
-
-#include "TCanvas.h"
-#include "TF1.h"
-#include "TH1.h"
-#include "TH1D.h"
-#include "TMath.h"
-#include "TStyle.h"
-
-using ROOT::RDataFrame;
-using namespace ROOT::VecOps;
-
-void emcal_barrel_electrons_analysis(const char* input_fname = "rec_emcal_barrel_uniform_electrons.root")
-{
-  // 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();
-  ROOT::RDataFrame d0("events", input_fname);
-
-  // Thrown Energy [GeV]
-  auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
-    std::vector<double> result;
-    result.push_back(TMath::Sqrt(input[2].psx * input[2].psx + input[2].psy * input[2].psy +
-                                 input[2].psz * input[2].psz + input[2].mass * input[2].mass));
-    return result;
-  };
-
-  // Reconstructed Energy [GeV] in XY merger
-  auto ErecXY = [](const std::vector<eic::CalorimeterHitData>& evt) {
-    std::vector<double> result;
-    auto                total_eng = 0.0;
-    for (const auto& i : evt)
-      total_eng += i.energy;
-    result.push_back(total_eng / 1.e+3);
-    return result;
-  };
-
-  // Reconstructed Energy [GeV] in Z merger
-  auto ErecZ = [](const std::vector<eic::CalorimeterHitData>& evt) {
-    std::vector<double> result;
-    auto                total_eng = 0.0;
-    for (const auto& i : evt)
-      total_eng += i.energy;
-    result.push_back(total_eng / 1.e+3);
-    return result;
-  };
-
-  // Number of Clusters
-  auto ncluster = [](const std::vector<eic::ClusterData>& evt) { return (int)evt.size(); };
-
-  // Cluster Energy [GeV]
-  auto Ecluster = [](const std::vector<eic::ClusterData>& evt) {
-    std::vector<double> result;
-    for (const auto& i : evt)
-      result.push_back(i.energy / 1.e+3);
-    return result;
-  };
-
-  // Sampling fraction = Esampling / Ethrown
-  auto fsam = [](const std::vector<double>& sampled, const std::vector<double>& thrown) {
-    std::vector<double> result;
-    for (const auto& E1 : thrown) {
-      for (const auto& E2 : sampled)
-        result.push_back(E2 / E1);
-    }
-    return result;
-  };
-
-  // Define variables
-  auto d1 = d0.Define("Ethr", Ethr, {"mcparticles2"})
-                .Define("ErecXY", ErecXY, {"RecoEcalBarrelHitsXY"})
-                .Define("ErecZ", ErecZ, {"RecoEcalBarrelHitsZ"})
-                .Define("ncluster", ncluster, {"EcalBarrelClusters"})
-                .Define("Ecluster", Ecluster, {"EcalBarrelClusters"})
-                .Define("fsam", fsam, {"Ecluster", "Ethr"});
-
-  // Define Histograms
-  auto hEthr   = d1.Histo1D({"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 6.5}, "Ethr");
-  auto hErecXY = d1.Histo1D(
-      {"hErecXY", "Reconstructed Energy in XY merger; Reconstructed Energy [GeV]; Events", 100, 0.0, 6.5}, "ErecXY");
-  auto hErecZ = d1.Histo1D(
-      {"hErecZ", "Reconstructed Energy in Z merger; Reconstructed Energy [GeV]; Events", 100, 0.0, 6.5}, "ErecZ");
-  auto hNCluster = d1.Histo1D({"hNCluster", "Number of Clusters; # of Clusters; Events", 100, 0.0, 100.0}, "ncluster");
-  auto hEcluster = d1.Histo1D({"hEcluster", "Cluster Energy; Cluster Energy [GeV]; Events", 100, 0.0, 6.5}, "Ecluster");
-  auto hfsam     = d1.Histo1D({"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1}, "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);
-  hEthr->GetYaxis()->SetTitleOffset(1.4);
-  hEthr->SetLineWidth(2);
-  hEthr->SetLineColor(kBlue);
-  hEthr->DrawClone();
-  c1->SaveAs("results/emcal_barrel_electrons_Ethr.png");
-  c1->SaveAs("results/emcal_barrel_electrons_Ethr.pdf");
-
-  TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
-  c2->SetLogy(1);
-  hErecXY->GetYaxis()->SetTitleOffset(1.4);
-  hErecXY->SetLineWidth(2);
-  hErecXY->SetLineColor(kBlue);
-  hErecXY->DrawClone();
-  c2->SaveAs("results/emcal_barrel_electrons_ErecXY.png");
-  c2->SaveAs("results/emcal_barrel_electrons_ErecXY.pdf");
-
-  TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
-  c3->SetLogy(1);
-  hErecZ->GetYaxis()->SetTitleOffset(1.4);
-  hErecZ->SetLineWidth(2);
-  hErecZ->SetLineColor(kBlue);
-  hErecZ->DrawClone();
-  c3->SaveAs("results/emcal_electrons_ErecZ.png");
-  c3->SaveAs("results/emcal_electrons_ErecZ.pdf");
-
-  TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
-  c4->SetLogy(1);
-  hNCluster->GetYaxis()->SetTitleOffset(1.6);
-  hNCluster->SetLineWidth(2);
-  hNCluster->SetLineColor(kBlue);
-  hNCluster->DrawClone();
-  c4->SaveAs("results/emcal_barrel_electrons_ncluster.png");
-  c4->SaveAs("results/emcal_barrel_electrons_ncluster.pdf");
-
-  TCanvas* c5 = new TCanvas("c5", "c5", 700, 500);
-  c5->SetLogy(1);
-  hEcluster->GetYaxis()->SetTitleOffset(1.4);
-  hEcluster->SetLineWidth(2);
-  hEcluster->SetLineColor(kBlue);
-  hEcluster->DrawClone();
-  c5->SaveAs("results/emcal_barrel_electrons_Ecluster.png");
-  c5->SaveAs("results/emcal_barrel_electrons_Ecluster.pdf");
-
-  TCanvas* c6 = new TCanvas("c6", "c6", 700, 500);
-  c6->SetLogy(1);
-  hfsam->GetYaxis()->SetTitleOffset(1.4);
-  hfsam->SetLineWidth(2);
-  hfsam->SetLineColor(kBlue);
-  hfsam->Fit("gaus", "", "", 0.01, 0.1);
-  // From emcal_barrel_pions_analysis.cxx
-  // Commented out for now as giving issues with new container (S. Joosten)
-  // hfsam->GetFunction("gaus")->SetLineWidth(2);
-  // hfsam->GetFunction("gaus")->SetLineColor(kRed);
-  hfsam->DrawClone();
-  c6->SaveAs("results/emcal_barrel_electrons_fsam.png");
-  c6->SaveAs("results/emcal_barrel_electrons_fsam.pdf");
-}
diff --git a/benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons_reader.cxx b/benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons_reader.cxx
deleted file mode 100644
index 75c2dae9..00000000
--- a/benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons_reader.cxx
+++ /dev/null
@@ -1,125 +0,0 @@
-//////////////////////////
-// EMCAL Barrel detector
-// Electron dataset
-// J.KIM 04/02/2021
-//////////////////////////
-#include "HepMC3/GenEvent.h"
-#include "HepMC3/Print.h"
-#include "HepMC3/ReaderAscii.h"
-#include "HepMC3/WriterAscii.h"
-
-#include "TH1F.h"
-#include "TStyle.h"
-#include <iostream>
-
-using namespace HepMC3;
-
-void emcal_barrel_electrons_reader(double e_start = 0.0, double e_end = 30.0, const char* in_fname = "./data/emcal_barrel_electrons.hepmc") {
-  // 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);
-
-  ReaderAscii hepmc_input(in_fname);
-  int events_parsed = 0;
-  GenEvent evt(Units::GEV, Units::MM);
-
-  // Histograms
-  TH1F* h_electrons_energy = new TH1F("h_electron_energy", "electron energy;E [GeV];Events",         100, -0.5, 30.5);
-  TH1F* h_electrons_eta    = new TH1F("h_electron_eta",    "electron #eta;#eta;Events",              100, -10.0, 10.0);
-  TH1F* h_electrons_theta  = new TH1F("h_electron_theta",  "electron #theta;#theta [degree];Events", 100, -0.5, 180.5);
-  TH1F* h_electrons_phi    = new TH1F("h_electron_phi",    "electron #phi;#phi [degree];Events",     100, -180.5, 180.5);
-  TH2F* h_electrons_pzpt   = new TH2F("h_electrons_pzpt",  "electron pt vs pz;pt [GeV];pz [GeV]",    100, -0.5, 30.5, 100, -30.5, 30.5);
-  TH2F* h_electrons_pxpy   = new TH2F("h_electrons_pxpy",  "electron px vs py;px [GeV];py [GeV]",    100, -30.5, 30.5, 100, -30.5, 30.5);
-  TH3F* h_electrons_p      = new TH3F("h_electron_p",      "electron p;px [GeV];py [GeV];pz [GeV]",  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;
-
-    for (const auto& v : evt.vertices()) {
-      for (const auto& p : v->particles_out()) {
-        if (p->pid() == 11) {
-          h_electrons_energy->Fill(p->momentum().e());
-          h_electrons_eta->Fill(p->momentum().eta());
-          h_electrons_theta->Fill(p->momentum().theta() * TMath::RadToDeg());
-          h_electrons_phi->Fill(p->momentum().phi() * TMath::RadToDeg());
-          h_electrons_pzpt->Fill(TMath::Sqrt(p->momentum().px() * p->momentum().px() + p->momentum().py() * p->momentum().py()), p->momentum().pz());
-          h_electrons_pxpy->Fill(p->momentum().px(), p->momentum().py());
-          h_electrons_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_electrons_energy->GetYaxis()->SetTitleOffset(1.8);
-  h_electrons_energy->SetLineWidth(2);
-  h_electrons_energy->SetLineColor(kBlue);
-  h_electrons_energy->DrawClone();
-  c->SaveAs("results/input_emcal_barrel_electrons_energy.png");
-  c->SaveAs("results/input_emcal_barrel_electrons_energy.pdf");
-
-  TCanvas* c1 = new TCanvas("c1", "c1", 500, 500);
-  h_electrons_eta->GetYaxis()->SetTitleOffset(1.9);
-  h_electrons_eta->SetLineWidth(2);
-  h_electrons_eta->SetLineColor(kBlue);
-  h_electrons_eta->DrawClone();
-  c1->SaveAs("results/input_emcal_barrel_electrons_eta.png");
-  c1->SaveAs("results/input_emcal_barrel_electrons_eta.pdf");
-
-  TCanvas* c2 = new TCanvas("c2", "c2", 500, 500);
-  h_electrons_theta->GetYaxis()->SetTitleOffset(1.8);
-  h_electrons_theta->SetLineWidth(2);
-  h_electrons_theta->SetLineColor(kBlue);
-  h_electrons_theta->DrawClone();
-  c2->SaveAs("results/input_emcal_barrel_electrons_theta.png");
-  c2->SaveAs("results/input_emcal_barrel_electrons_theta.pdf");
-
-  TCanvas* c3 = new TCanvas("c3", "c3", 500, 500);
-  h_electrons_phi->GetYaxis()->SetTitleOffset(1.8);
-  h_electrons_phi->SetLineWidth(2);
-  h_electrons_phi->GetYaxis()->SetRangeUser(0.0, h_electrons_phi->GetMaximum() + 100.0);
-  h_electrons_phi->SetLineColor(kBlue);
-  h_electrons_phi->DrawClone();
-  c3->SaveAs("results/input_emcal_barrel_electrons_phi.png");
-  c3->SaveAs("results/input_emcal_barrel_electrons_phi.pdf");
-
-  TCanvas* c4 = new TCanvas("c4", "c4", 500, 500);
-  h_electrons_pzpt->GetYaxis()->SetTitleOffset(1.4);
-  h_electrons_pzpt->SetLineWidth(2);
-  h_electrons_pzpt->SetLineColor(kBlue);
-  h_electrons_pzpt->DrawClone("COLZ");
-  c4->SaveAs("results/input_emcal_barrel_electrons_pzpt.png");
-  c4->SaveAs("results/input_emcal_barrel_electrons_pzpt.pdf");
-
-  TCanvas* c5 = new TCanvas("c5", "c5", 500, 500);
-  h_electrons_pxpy->GetYaxis()->SetTitleOffset(1.4);
-  h_electrons_pxpy->SetLineWidth(2);
-  h_electrons_pxpy->SetLineColor(kBlue);
-  h_electrons_pxpy->DrawClone("COLZ");
-  c5->SaveAs("results/input_emcal_barrel_electrons_pxpy.png");
-  c5->SaveAs("results/input_emcal_barrel_electrons_pxpy.pdf");
-
-  TCanvas* c6 = new TCanvas("c6", "c6", 500, 500);
-  h_electrons_p->GetYaxis()->SetTitleOffset(1.8);
-  h_electrons_p->GetXaxis()->SetTitleOffset(1.6);
-  h_electrons_p->GetZaxis()->SetTitleOffset(1.6);
-  h_electrons_p->SetLineWidth(2);
-  h_electrons_p->SetLineColor(kBlue);
-  h_electrons_p->DrawClone();
-  c6->SaveAs("results/input_emcal_barrel_electrons_p.png");
-  c6->SaveAs("results/input_emcal_barrel_electrons_p.pdf");
-}
diff --git a/benchmarks/ecal/scripts/deprecated/emcal_barrel_pions.cxx b/benchmarks/ecal/scripts/deprecated/emcal_barrel_pions.cxx
deleted file mode 100644
index 02537b56..00000000
--- a/benchmarks/ecal/scripts/deprecated/emcal_barrel_pions.cxx
+++ /dev/null
@@ -1,80 +0,0 @@
-//////////////////////////////////////////////////////////////
-// EMCAL Barrel detector
-// Single Pion dataset
-// J.KIM 04/04/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>
-
-using namespace HepMC3;
-
-void emcal_barrel_pions(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/emcal_barrel_pions.hepmc") {
-  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 111 - pi0
-    // 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);
-
-    // Generates random vectors, uniformly distributed over the surface of a
-    // sphere of given radius, in this case momentum.
-    // r1->Sphere(px, py, pz, p);
-
-    // type 1 is final state
-    // pdgid 211 - pion+ 139.570 MeV/c^2
-    GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.139570 * 0.139570))), 211, 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/ecal/scripts/deprecated/emcal_barrel_pions_analysis.cxx b/benchmarks/ecal/scripts/deprecated/emcal_barrel_pions_analysis.cxx
deleted file mode 100644
index b857a388..00000000
--- a/benchmarks/ecal/scripts/deprecated/emcal_barrel_pions_analysis.cxx
+++ /dev/null
@@ -1,174 +0,0 @@
-////////////////////////////////////////
-// Read reconstruction ROOT output file
-// Plot variables
-////////////////////////////////////////
-
-#include "ROOT/RDataFrame.hxx"
-#include <iostream>
-
-#include "dd4pod/CalorimeterHitCollection.h"
-#include "dd4pod/Geant4ParticleCollection.h"
-#include "dd4pod/TrackerHitCollection.h"
-#include "eicd/CalorimeterHitCollection.h"
-#include "eicd/CalorimeterHitData.h"
-#include "eicd/ClusterCollection.h"
-#include "eicd/ClusterData.h"
-#include "eicd/RawCalorimeterHitCollection.h"
-#include "eicd/RawCalorimeterHitData.h"
-
-#include "TCanvas.h"
-#include "TF1.h"
-#include "TH1.h"
-#include "TH1D.h"
-#include "TMath.h"
-#include "TStyle.h"
-
-using ROOT::RDataFrame;
-using namespace ROOT::VecOps;
-
-void emcal_barrel_pions_analysis(const char* input_fname = "rec_emcal_barrel_uniform_pions.root")
-{
-  // 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();
-  ROOT::RDataFrame d0("events", input_fname);
-
-  // Thrown Energy [GeV]
-  auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
-    std::vector<double> result;
-    result.push_back(TMath::Sqrt(input[2].psx * input[2].psx + input[2].psy * input[2].psy +
-                                 input[2].psz * input[2].psz + input[2].mass * input[2].mass));
-    return result;
-  };
-
-  // Reconstructed Energy [GeV] in XY merger
-  auto ErecXY = [](const std::vector<eic::CalorimeterHitData>& evt) {
-    std::vector<double> result;
-    auto                total_eng = 0.0;
-    for (const auto& i : evt)
-      total_eng += i.energy;
-    result.push_back(total_eng / 1.e+3);
-    return result;
-  };
-
-  // Reconstructed Energy [GeV] in Z merger
-  auto ErecZ = [](const std::vector<eic::CalorimeterHitData>& evt) {
-    std::vector<double> result;
-    auto                total_eng = 0.0;
-    for (const auto& i : evt)
-      total_eng += i.energy;
-    result.push_back(total_eng / 1.e+3);
-    return result;
-  };
-
-  // Number of Clusters
-  auto ncluster = [](const std::vector<eic::ClusterData>& evt) { return (int)evt.size(); };
-
-  // Cluster Energy [GeV]
-  auto Ecluster = [](const std::vector<eic::ClusterData>& evt) {
-    std::vector<double> result;
-    for (const auto& i : evt)
-      result.push_back(i.energy / 1.e+3);
-    return result;
-  };
-
-  // Sampling fraction = Esampling / Ethrown
-  auto fsam = [](const std::vector<double>& sampled, const std::vector<double>& thrown) {
-    std::vector<double> result;
-    for (const auto& E1 : thrown) {
-      for (const auto& E2 : sampled)
-        result.push_back(E2 / E1);
-    }
-    return result;
-  };
-
-  // Define variables
-  auto d1 = d0.Define("Ethr", Ethr, {"mcparticles2"})
-                .Define("ErecXY", ErecXY, {"RecoEcalBarrelHitsXY"})
-                .Define("ErecZ", ErecZ, {"RecoEcalBarrelHitsZ"})
-                .Define("ncluster", ncluster, {"EcalBarrelClusters"})
-                .Define("Ecluster", Ecluster, {"EcalBarrelClusters"})
-                .Define("fsam", fsam, {"Ecluster", "Ethr"});
-
-  // Define Histograms
-  auto hEthr   = d1.Histo1D({"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 6.5}, "Ethr");
-  auto hErecXY = d1.Histo1D(
-      {"hErecXY", "Reconstructed Energy in XY merger; Reconstructed Energy [GeV]; Events", 100, 0.0, 6.5}, "ErecXY");
-  auto hErecZ = d1.Histo1D(
-      {"hErecZ", "Reconstructed Energy in Z merger; Reconstructed Energy [GeV]; Events", 100, 0.0, 6.5}, "ErecZ");
-  auto hNCluster = d1.Histo1D({"hNCluster", "Number of Clusters; # of Clusters; Events", 100, 0.0, 100.0}, "ncluster");
-  auto hEcluster = d1.Histo1D({"hEcluster", "Cluster Energy; Cluster Energy [GeV]; Events", 100, 0.0, 6.5}, "Ecluster");
-  auto hfsam     = d1.Histo1D({"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1}, "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);
-  hEthr->GetYaxis()->SetTitleOffset(1.4);
-  hEthr->SetLineWidth(2);
-  hEthr->SetLineColor(kBlue);
-  hEthr->DrawClone();
-  c1->SaveAs("results/emcal_barrel_pions_Ethr.png");
-  c1->SaveAs("results/emcal_barrel_pions_Ethr.pdf");
-
-  TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
-  c2->SetLogy(1);
-  hErecXY->GetYaxis()->SetTitleOffset(1.4);
-  hErecXY->SetLineWidth(2);
-  hErecXY->SetLineColor(kBlue);
-  hErecXY->DrawClone();
-  c2->SaveAs("results/emcal_barrel_pions_ErecXY.png");
-  c2->SaveAs("results/emcal_barrel_pions_ErecXY.pdf");
-
-  TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
-  c3->SetLogy(1);
-  hErecZ->GetYaxis()->SetTitleOffset(1.4);
-  hErecZ->SetLineWidth(2);
-  hErecZ->SetLineColor(kBlue);
-  hErecZ->DrawClone();
-  c3->SaveAs("results/emcal_pions_ErecZ.png");
-  c3->SaveAs("results/emcal_pions_ErecZ.pdf");
-
-  TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
-  c4->SetLogy(1);
-  hNCluster->GetYaxis()->SetTitleOffset(1.6);
-  hNCluster->SetLineWidth(2);
-  hNCluster->SetLineColor(kBlue);
-  hNCluster->DrawClone();
-  c4->SaveAs("results/emcal_barrel_pions_ncluster.png");
-  c4->SaveAs("results/emcal_barrel_pions_ncluster.pdf");
-
-  TCanvas* c5 = new TCanvas("c5", "c5", 700, 500);
-  c5->SetLogy(1);
-  hEcluster->GetYaxis()->SetTitleOffset(1.4);
-  hEcluster->SetLineWidth(2);
-  hEcluster->SetLineColor(kBlue);
-  hEcluster->DrawClone();
-  c5->SaveAs("results/emcal_barrel_pions_Ecluster.png");
-  c5->SaveAs("results/emcal_barrel_pions_Ecluster.pdf");
-
-  TCanvas* c6 = new TCanvas("c6", "c6", 700, 500);
-  c6->SetLogy(1);
-  hfsam->GetYaxis()->SetTitleOffset(1.4);
-  hfsam->SetLineWidth(2);
-  hfsam->SetLineColor(kBlue);
-  hfsam->Fit("gaus", "", "", 0.01, 0.1);
-  // Commented out for now as giving issues with new container (S. Joosten)
-  // hfsam->GetFunction("gaus")->SetLineWidth(2);
-  // hfsam->GetFunction("gaus")->SetLineColor(kRed);
-  hfsam->DrawClone();
-  c6->SaveAs("results/emcal_barrel_pions_fsam.png");
-  c6->SaveAs("results/emcal_barrel_pions_fsam.pdf");
-}
diff --git a/benchmarks/ecal/scripts/deprecated/emcal_barrel_pions_reader.cxx b/benchmarks/ecal/scripts/deprecated/emcal_barrel_pions_reader.cxx
deleted file mode 100644
index e0377148..00000000
--- a/benchmarks/ecal/scripts/deprecated/emcal_barrel_pions_reader.cxx
+++ /dev/null
@@ -1,125 +0,0 @@
-//////////////////////////
-// EMCAL Barrel detector
-// Pion dataset
-// J.KIM 04/04/2021
-//////////////////////////
-#include "HepMC3/GenEvent.h"
-#include "HepMC3/Print.h"
-#include "HepMC3/ReaderAscii.h"
-#include "HepMC3/WriterAscii.h"
-
-#include "TH1F.h"
-#include "TStyle.h"
-#include <iostream>
-
-using namespace HepMC3;
-
-void emcal_barrel_pions_reader(double e_start = 0.0, double e_end = 30.0, const char* in_fname = "./data/emcal_barrel_pions.hepmc") {
-  // 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);
-
-  ReaderAscii hepmc_input(in_fname);
-  int events_parsed = 0;
-  GenEvent evt(Units::GEV, Units::MM);
-
-  // Histograms
-  TH1F* h_pions_energy = new TH1F("h_pions_energy", "pion energy;E [GeV];Events",         100, -0.5, 30.5);
-  TH1F* h_pions_eta    = new TH1F("h_pions_eta",    "pion #eta;#eta;Events",              100, -10.0, 10.0);
-  TH1F* h_pions_theta  = new TH1F("h_pions_theta",  "pion #theta;#theta [degree];Events", 100, -0.5, 180.5);
-  TH1F* h_pions_phi    = new TH1F("h_pions_phi",    "pion #phi;#phi [degree];Events",     100, -180.5, 180.5);
-  TH2F* h_pions_pzpt   = new TH2F("h_pions_pzpt",   "pion pt vs pz;pt [GeV];pz [GeV]",    100, -0.5, 30.5, 100, -30.5, 30.5);
-  TH2F* h_pions_pxpy   = new TH2F("h_pions_pxpy",   "pion px vs py;px [GeV];py [GeV]",    100, -30.5, 30.5, 100, -30.5, 30.5);
-  TH3F* h_pions_p      = new TH3F("h_pions_p",      "pion p;px [GeV];py [GeV];pz [GeV]",  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;
-
-    for (const auto& v : evt.vertices()) {
-      for (const auto& p : v->particles_out()) {
-        if (p->pid() == 11) {
-          h_pions_energy->Fill(p->momentum().e());
-          h_pions_eta->Fill(p->momentum().eta());
-          h_pions_theta->Fill(p->momentum().theta() * TMath::RadToDeg());
-          h_pions_phi->Fill(p->momentum().phi() * TMath::RadToDeg());
-          h_pions_pzpt->Fill(TMath::Sqrt(p->momentum().px() * p->momentum().px() + p->momentum().py() * p->momentum().py()), p->momentum().pz());
-          h_pions_pxpy->Fill(p->momentum().px(), p->momentum().py());
-          h_pions_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_pions_energy->GetYaxis()->SetTitleOffset(1.8);
-  h_pions_energy->SetLineWidth(2);
-  h_pions_energy->SetLineColor(kBlue);
-  h_pions_energy->DrawClone();
-  c->SaveAs("results/input_emcal_barrel_pions_energy.png");
-  c->SaveAs("results/input_emcal_barrel_pions_energy.pdf");
-
-  TCanvas* c1 = new TCanvas("c1", "c1", 500, 500);
-  h_pions_eta->GetYaxis()->SetTitleOffset(1.9);
-  h_pions_eta->SetLineWidth(2);
-  h_pions_eta->SetLineColor(kBlue);
-  h_pions_eta->DrawClone();
-  c1->SaveAs("results/input_emcal_barrel_pions_eta.png");
-  c1->SaveAs("results/input_emcal_barrel_pions_eta.pdf");
-
-  TCanvas* c2 = new TCanvas("c2", "c2", 500, 500);
-  h_pions_theta->GetYaxis()->SetTitleOffset(1.8);
-  h_pions_theta->SetLineWidth(2);
-  h_pions_theta->SetLineColor(kBlue);
-  h_pions_theta->DrawClone();
-  c2->SaveAs("results/input_emcal_barrel_pions_theta.png");
-  c2->SaveAs("results/input_emcal_barrel_pions_theta.pdf");
-
-  TCanvas* c3 = new TCanvas("c3", "c3", 500, 500);
-  h_pions_phi->GetYaxis()->SetTitleOffset(1.8);
-  h_pions_phi->SetLineWidth(2);
-  h_pions_phi->GetYaxis()->SetRangeUser(0.0, h_pions_phi->GetMaximum() + 100.0);
-  h_pions_phi->SetLineColor(kBlue);
-  h_pions_phi->DrawClone();
-  c3->SaveAs("results/input_emcal_barrel_pions_phi.png");
-  c3->SaveAs("results/input_emcal_barrel_pions_phi.pdf");
-
-  TCanvas* c4 = new TCanvas("c4", "c4", 500, 500);
-  h_pions_pzpt->GetYaxis()->SetTitleOffset(1.4);
-  h_pions_pzpt->SetLineWidth(2);
-  h_pions_pzpt->SetLineColor(kBlue);
-  h_pions_pzpt->DrawClone("COLZ");
-  c4->SaveAs("results/input_emcal_barrel_pions_pzpt.png");
-  c4->SaveAs("results/input_emcal_barrel_pions_pzpt.pdf");
-
-  TCanvas* c5 = new TCanvas("c5", "c5", 500, 500);
-  h_pions_pxpy->GetYaxis()->SetTitleOffset(1.4);
-  h_pions_pxpy->SetLineWidth(2);
-  h_pions_pxpy->SetLineColor(kBlue);
-  h_pions_pxpy->DrawClone("COLZ");
-  c5->SaveAs("results/input_emcal_barrel_pions_pxpy.png");
-  c5->SaveAs("results/input_emcal_barrel_pions_pxpy.pdf");
-
-  TCanvas* c6 = new TCanvas("c6", "c6", 500, 500);
-  h_pions_p->GetYaxis()->SetTitleOffset(1.8);
-  h_pions_p->GetXaxis()->SetTitleOffset(1.6);
-  h_pions_p->GetZaxis()->SetTitleOffset(1.6);
-  h_pions_p->SetLineWidth(2);
-  h_pions_p->SetLineColor(kBlue);
-  h_pions_p->DrawClone();
-  c6->SaveAs("results/input_emcal_barrel_pions_p.png");
-  c6->SaveAs("results/input_emcal_barrel_pions_p.pdf");
-}
diff --git a/benchmarks/ecal/scripts/deprecated/emcal_electrons.cxx b/benchmarks/ecal/scripts/deprecated/emcal_electrons.cxx
deleted file mode 100644
index 82296521..00000000
--- a/benchmarks/ecal/scripts/deprecated/emcal_electrons.cxx
+++ /dev/null
@@ -1,90 +0,0 @@
-//////////////////////////////////////////////////////////////
-// Crystal EMCAL detector 
-// Single Electron dataset
-// J.KIM 07/20/2020
-//
-// 10/4/2020 
-// Changed to have isotropic distribution in momentum sphere
-//////////////////////////////////////////////////////////////
-#include "HepMC3/GenEvent.h"
-#include "HepMC3/ReaderAscii.h"
-#include "HepMC3/WriterAscii.h"
-#include "HepMC3/Print.h"
-
-#include <iostream>
-#include<random>
-#include<cmath>
-#include <math.h>
-#include <TMath.h>
-
-using namespace HepMC3;
-
-void emcal_electrons(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/emcal_electron_0GeVto30GeV_1000kEvt.hepmc")
-{
-  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 * (160.0 / 180.0));
-  double cos_theta_max = std::cos(M_PI * (175.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 111 - pi0
-    // 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
-    // Momentum starting at 0 GeV/c to 30 GeV/c
-    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);
-
-    // Generates random vectors, uniformly distributed over the surface of a
-    // sphere of given radius, in this case momentum.
-    //r1->Sphere(px, py, pz, p);
-
-    // type 1 is final state
-    // pdgid 11 - electron 0.510 MeV/c^2
-    GenParticlePtr p3 = std::make_shared<GenParticle>(
-        FourVector(
-            px, py, pz,
-            sqrt(p*p + (0.000511 * 0.000511))),
-        11, 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/ecal/scripts/deprecated/emcal_electrons_analysis.cxx b/benchmarks/ecal/scripts/deprecated/emcal_electrons_analysis.cxx
deleted file mode 100644
index 1d8a20c6..00000000
--- a/benchmarks/ecal/scripts/deprecated/emcal_electrons_analysis.cxx
+++ /dev/null
@@ -1,492 +0,0 @@
-////////////////////////////////////////
-// Read reconstruction ROOT output file
-// Plot variables
-////////////////////////////////////////
-
-#include "ROOT/RDataFrame.hxx"
-#include <iostream>
-
-#include "dd4pod/Geant4ParticleCollection.h"
-#include "eicd/ClusterCollection.h"
-#include "eicd/ClusterData.h"
-
-#include "TCanvas.h"
-#include "TStyle.h"
-#include "TMath.h"
-#include "TH1.h"
-#include "TH1D.h"
-
-// If root cannot find a dd4pod header make sure to add the include dir to ROOT_INCLUDE_PATH:
-//  export ROOT_INCLUDE_PATH=$HOME/include:$ROOT_INCLUDE_PATH
-//  export ROOT_INCLUDE_PATH=/home/jihee/stow/development/include:$ROOT_INCLUDE_PATH
-
-using ROOT::RDataFrame;
-using namespace ROOT::VecOps;
-
-void emcal_electrons_analysis(const char* input_fname = "rec_emcal_uniform_electrons.root")
-{
-  // 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();
-  ROOT::RDataFrame d0("events", input_fname);
-
-  // Number of Clusters
-  auto ncluster = [] (const std::vector<eic::ClusterData>& evt) {return (int) evt.size(); };
-
-  // Cluster Energy [GeV]
-  auto clusterE = [] (const std::vector<eic::ClusterData>& evt) {
-  	std::vector<float> result;
-  	for (const auto& i: evt)
-    		result.push_back(i.energy);
-  return result;
-  };
-
-  // Scattering Angle [degree]
-  auto theta = [] (const std::vector<eic::ClusterData>& evt) {
-	std::vector<float> result;
-	for(const auto& i: evt) {
-		auto r = TMath::Sqrt(i.position.x*i.position.x + i.position.y*i.position.y + i.position.z*i.position.z);
-		result.push_back(TMath::ACos(i.position.z/r)*TMath::RadToDeg());
-	}
-  return result;
-  };
-
-  // Pseudo-rapidity
-  auto eta = [] (const std::vector<eic::ClusterData>& evt) {
-	std::vector<float> result;
-	for(const auto& i: evt) {
-		auto r  = TMath::Sqrt(i.position.x*i.position.x + i.position.y*i.position.y + i.position.z*i.position.z);
-		auto th = TMath::ACos(i.position.z/r)*TMath::RadToDeg(); 
-		result.push_back(-1.0 * TMath::Log(TMath::Tan((th*TMath::DegToRad())/2.0)));
-	}
-  return result;
-  };
-
-  // Mass [GeV]
-  auto mass2 = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
-  	std::vector<float> result;
-    	result.push_back(input[2].mass*input[2].mass);
-  return result;
-  };
-
-  // Momentum [GeV]
-  auto p_rec = [](const std::vector<float>& energy_term, const std::vector<float>& mass_term) {
-	std::vector<float> result;
-	for (const auto& e : energy_term) {
-		for (const auto& m2 : mass_term) {
-			result.push_back(TMath::Sqrt((e*e) - m2));
-		}
-	}
-  return result;
-  };
-  
-  // Thrown Momentum [GeV]
-  auto p_thr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
-        std::vector<float> result;
-        result.push_back(TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz));
-  return result;
-  };
-
-  // Thrown Energy [GeV]
-  auto E_thr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
-  	std::vector<float> result;
-    	result.push_back(TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass));
-  return result;
-  };
-
-  // Energy Resolution
-  auto E_res = [] (const std::vector<float>& Erec, const std::vector<float>& Ethr) {
-	std::vector<float> result;
-  	for (const auto& E1 : Ethr) {
-    		for (const auto& E2 : Erec) {
-      			result.push_back((E2-E1)/E1);
-    		}
-  	}
-  return result;
-  };
-
-  // Momentum Ratio
-  auto p_ratio = [] (const std::vector<float>& Prec, const std::vector<float>& Pthr) {
-	std::vector<float> result;
-  	for (const auto& P1 : Pthr) {
-    		for (const auto& P2 : Prec) {
-      			result.push_back(P2/P1);
-    		}
-  	}
-  return result;
-  };
-  
-  // Define variables
-  auto d1 = d0.Define("ncluster", ncluster, {"EcalClusters"})
-	      .Define("clusterE", clusterE, {"EcalClusters"})
-	      .Define("theta",    theta,    {"EcalClusters"})
-	      .Define("eta",      eta,      {"EcalClusters"})
-	      .Define("mass2",    mass2,    {"mcparticles2"})
-	      .Define("p_rec",    p_rec,    {"clusterE","mass2"})
-	      .Define("p_thr",    p_thr,    {"mcparticles2"})
-	      .Define("E_thr",    E_thr,    {"mcparticles2"})
-	      .Define("E_res",    E_res,    {"clusterE","E_thr"})
-	      .Define("p_ratio",  p_ratio,  {"p_rec","p_thr"})
-	      ;
-
-  // Define Histograms
-  auto hNCluster = d1.Histo1D({"hNCluster", "Number of Clusters; # of Clusters; Events",    5, -0.5, 4.5},      "ncluster");
-  auto hClusterE = d1.Histo1D({"hClusterE", "Cluster Energy; Cluster Energy [GeV]; Events", 100, -0.5, 30.5},   "clusterE");
-  auto hTheta    = d1.Histo1D({"hTheta",    "Scattering Angle; #theta [degree]; Events",    100, 130.0, 180.0}, "theta");
-  auto hEta      = d1.Histo1D({"hEta",      "Pseudo-rapidity; #eta; Events",                100, -5.0, 0.0},    "eta");
-  auto hPrec     = d1.Histo1D({"hPrec",     "Reconstructed Momentum; p_{rec}[GeV]; Events", 100, -0.5, 30.5},   "p_rec");
-  auto hPthr     = d1.Histo1D({"hPthr",     "Thrown Momentum; p_{thr}[GeV]; Events",        100, -0.5, 30.5},   "p_thr");
-  auto hEthr     = d1.Histo1D({"hEthr",     "Thrown Energy; E_{thr}[GeV]; Events",          100, -0.5, 30.5},   "E_thr");
-
-  // Select Events with One Cluster
-  auto d2 = d1.Filter("ncluster==1");
-  auto hClusterE1       = d2.Histo1D({"hClusterE1", "One Cluster Energy; Cluster Energy [GeV]; Events", 100, -0.5, 30.5}, "clusterE");
-  auto hEres            = d2.Histo1D({"hEres",      "Energy Resolution; #DeltaE/E; Events",             100, -1.0, 1.0},  "E_res");
-  auto hEres_lmom       = d2.Filter([=](std::vector<dd4pod::Geant4ParticleData> const& input) {
-                                        auto Pthr = TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass);
-                                        if (Pthr <= 15.0)
-                                                return true;
-                                        return false;}, {"mcparticles2"})
-                            .Histo1D({"hEres_lmom", "Energy Resolution; #DeltaE/E; Events",             100, -1.0, 1.0},  "E_res");
-  auto hEres_hmom       = d2.Filter([=](std::vector<dd4pod::Geant4ParticleData> const& input) {
-                                        auto Pthr = TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass);
-                                        if (Pthr > 15.0)
-                                                return true;
-                                        return false;}, {"mcparticles2"})
-                            .Histo1D({"hEres_hmom", "Energy Resolution; #DeltaE/E; Events",             100, -1.0, 1.0},  "E_res");
-  auto hEres_leta       = d2.Filter([=](const std::vector<eic::ClusterData>& evt) {
-                                        for(const auto& i: evt) {
-		  				auto r   = TMath::Sqrt(i.position.x*i.position.x + i.position.y*i.position.y + i.position.z*i.position.z);
-						auto th  = TMath::ACos(i.position.z/r)*TMath::RadToDeg();
-						auto eta = -1.0 * TMath::Log(TMath::Tan((th*TMath::DegToRad())/2.0));
-						if (eta <= -2.25)
-                                                	return true;
-					}
-                                        return false;}, {"EcalClusters"})
-                            .Histo1D({"hEres_leta", "Energy Resolution; #DeltaE/E; Events",             100, -1.0, 1.0},  "E_res");
-  auto hEres_heta       = d2.Filter([=](const std::vector<eic::ClusterData>& evt) {
-					for(const auto& i: evt) {
-                                                auto r   = TMath::Sqrt(i.position.x*i.position.x + i.position.y*i.position.y + i.position.z*i.position.z);
-                                                auto th  = TMath::ACos(i.position.z/r)*TMath::RadToDeg();
-                                                auto eta = -1.0 * TMath::Log(TMath::Tan((th*TMath::DegToRad())/2.0));
-                                                if (eta > -2.25)
-                                                        return true;
-                                        }
-                                        return false;}, {"EcalClusters"})
-                            .Histo1D({"hEres_heta", "Energy Resolution; #DeltaE/E; Events",             100, -1.0, 1.0},  "E_res");
-  auto hPratio          = d2.Histo1D({"hPratio",    "Momentum ratio; p_{rec}/p_{thr}; Events",          100, 0.0, 1.0},   "p_ratio");
-  auto hPthr_accepted   = d2.Filter([=] (const std::vector<float>& Prec, const std::vector<float>& Pthr) {
-  					for (const auto& P1 : Pthr) {
-    						for (const auto& P2 : Prec) {
-      							auto Pratio = P2/P1;
-							if (Pratio > 0.70) {
-								return true;
-							}
-    						}
-  					}
-  			           return false;}, {"p_rec","p_thr"})
-	  		    .Histo1D({"hPthr_accepted", "Thrown momentum for reconstructed particle; p_{thr} [GeV]; Events", 100, -0.5, 30.5}, "p_thr");
-  // 2D histograms
-  auto hPrec_Pthr       = d2.Histo2D({"hPrec_Pthr", "p_{rec} vs p_{thr}; p_{rec}[GeV]; p_{thr[GeV]}",     100, -0.5, 30.5, 100, -0.5, 30.5}, "p_rec", "p_thr");
-  auto hEres_Ethr       = d2.Histo2D({"hEres_Ethr", "Energy Resolution vs Ethr; #DeltaE/E; E_{thr}[GeV]", 100, -1.0, 1.0, 100, -0.5, 30.5},  "E_res", "E_thr");
-
-  // Cut on Radial Distance
-  auto d3 = d2.Filter([=] (const std::vector<eic::ClusterData>& evt) {
-		  			for(const auto& i: evt) {
-						auto pos_x = i.position.x;
-						auto pos_y = i.position.y;
-						auto radial_dist = TMath::Sqrt(pos_x*pos_x + pos_y*pos_y);
-						if (radial_dist > 18.0 && radial_dist < 30.0)
-							return true;
-					}
-		  		   return false;}, {"EcalClusters"});
-  auto hEres_cut      = d3.Histo1D({"hEres_cut",  "Energy Resolution; #DeltaE/E; Events",      100, -1.0, 1.0},    "E_res");
-  auto hEres_cut_lmom = d3.Filter([=](std::vector<dd4pod::Geant4ParticleData> const& input) {
-                                        auto Pthr = TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass);
-                                        if (Pthr <= 15.0)
-                                                return true;
-                                        return false;}, {"mcparticles2"})
-                            .Histo1D({"hEres_cut_lmom", "Energy Resolution; #DeltaE/E; Events",             100, -1.0, 1.0},  "E_res");
-  auto hEres_cut_hmom = d3.Filter([=](std::vector<dd4pod::Geant4ParticleData> const& input) {
-                                        auto Pthr = TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass);
-                                        if (Pthr > 15.0)
-                                                return true;
-                                        return false;}, {"mcparticles2"})
-                            .Histo1D({"hEres_cut_hmom", "Energy Resolution; #DeltaE/E; Events",             100, -1.0, 1.0},  "E_res");
-  auto hEres_cut_leta = d3.Filter([=](const std::vector<eic::ClusterData>& evt) {
-                                        for(const auto& i: evt) {
-                                                auto r   = TMath::Sqrt(i.position.x*i.position.x + i.position.y*i.position.y + i.position.z*i.position.z);
-                                                auto th  = TMath::ACos(i.position.z/r)*TMath::RadToDeg();
-                                                auto eta = -1.0 * TMath::Log(TMath::Tan((th*TMath::DegToRad())/2.0));
-                                                if (eta <= -2.25)
-                                                        return true;
-                                        }
-                                        return false;}, {"EcalClusters"})
-                            .Histo1D({"hEres_cut_leta", "Energy Resolution; #DeltaE/E; Events",             100, -1.0, 1.0},  "E_res");
-  auto hEres_cut_heta = d3.Filter([=](const std::vector<eic::ClusterData>& evt) {
-                                        for(const auto& i: evt) {
-                                                auto r   = TMath::Sqrt(i.position.x*i.position.x + i.position.y*i.position.y + i.position.z*i.position.z);
-                                                auto th  = TMath::ACos(i.position.z/r)*TMath::RadToDeg();
-                                                auto eta = -1.0 * TMath::Log(TMath::Tan((th*TMath::DegToRad())/2.0));
-                                                if (eta > -2.25)
-                                                        return true;
-                                        }
-                                        return false;}, {"EcalClusters"})
-                            .Histo1D({"hEres_cut_heta", "Energy Resolution; #DeltaE/E; Events",             100, -1.0, 1.0},  "E_res");
-  auto hTheta_cut     = d3.Histo1D({"hTheta_cut", "Scattering Angle; #theta [degree]; Events", 100, 130.0, 180.0}, "theta");
-  auto hEta_cut       = d3.Histo1D({"hEta_cut",   "Pseudo-rapidity; #eta; Events",             100, -5.0, 0.0},    "eta");
-
-  // 2D histograms
-  auto hPrec_Pthr_cut = d3.Histo2D({"hPrec_Pthr_cut", "p_{rec} vs p_{thr}; p_{rec}[GeV]; p_{thr[GeV]}",     100, -0.5, 30.5, 100, -0.5, 30.5}, "p_rec", "p_thr");
-  auto hEres_Ethr_cut = d3.Histo2D({"hEres_Ethr_cut", "Energy Resolution vs Ethr; #DeltaE/E; E_{thr}[GeV]", 100, -1.0, 1.0, 100, -0.5, 30.5},  "E_res", "E_thr");
-
-  // Event Counts
-  auto nevents_thrown    = d1.Count();
-  auto nevents_cluster1  = d2.Count();
-
-  std::cout << "Number of Events: " << (*nevents_cluster1) << " / " << (*nevents_thrown) << " = single cluster events/all \n";
-
-  // Draw Histograms
-  TCanvas *c1 = new TCanvas("c1", "c1", 500, 500);
-  c1->SetLogy(1); 
-  hNCluster->GetYaxis()->SetTitleOffset(1.6);
-  hNCluster->SetLineWidth(2);
-  hNCluster->SetLineColor(kBlue);
-  hNCluster->DrawClone();
-  c1->SaveAs("results/emcal_electrons_ncluster.png");
-  c1->SaveAs("results/emcal_electrons_ncluster.pdf");
-
-  TCanvas *c2 = new TCanvas("c2", "c2", 500, 500); 
-  hClusterE->GetYaxis()->SetTitleOffset(1.4);
-  hClusterE->SetLineWidth(2);
-  hClusterE->SetLineColor(kBlue);
-  hClusterE->DrawClone();
-  c2->SaveAs("results/emcal_electrons_clusterE.png");
-  c2->SaveAs("results/emcal_electrons_clusterE.pdf");
-
-  TCanvas *c3 = new TCanvas("c3", "c3", 500, 500);
-  hTheta->GetYaxis()->SetTitleOffset(1.4);
-  hTheta->SetLineWidth(2);
-  hTheta->SetLineColor(kBlue);
-  hTheta->DrawClone();
-  c3->SaveAs("results/emcal_electrons_theta.png"); 
-  c3->SaveAs("results/emcal_electrons_theta.pdf");
-
-  TCanvas *c4 = new TCanvas("c4", "c4", 500, 500);
-  hEta->GetYaxis()->SetTitleOffset(1.4);
-  hEta->SetLineWidth(2);
-  hEta->SetLineColor(kBlue);
-  hEta->DrawClone();
-  c4->SaveAs("results/emcal_electrons_eta.png");
-  c4->SaveAs("results/emcal_electrons_eta.pdf");
-
-  TCanvas *c5 = new TCanvas("c5", "c5", 500, 500);
-  hPrec->GetYaxis()->SetTitleOffset(1.4);
-  hPrec->SetLineWidth(2);
-  hPrec->SetLineColor(kBlue);
-  hPrec->DrawClone();
-  c5->SaveAs("results/emcal_electrons_Prec.png");
-  c5->SaveAs("results/emcal_electrons_Prec.pdf");
-
-  TCanvas *c6 = new TCanvas("c6", "c6", 500, 500);
-  hPthr->GetYaxis()->SetTitleOffset(1.4);
-  hPthr->SetLineWidth(2);
-  hPthr->SetLineColor(kBlue);
-  hPthr->DrawClone();
-  c6->SaveAs("results/emcal_electrons_Pthr.png");
-  c6->SaveAs("results/emcal_electrons_Pthr.pdf");
-
-  TCanvas *c7 = new TCanvas("c7", "c7", 500, 500);
-  hEthr->GetYaxis()->SetTitleOffset(1.4);
-  hEthr->SetLineWidth(2);
-  hEthr->SetLineColor(kBlue);
-  hEthr->DrawClone();
-  c7->SaveAs("results/emcal_electrons_Ethr.png");
-  c7->SaveAs("results/emcal_electrons_Ethr.pdf");
-
-  TCanvas *c8 = new TCanvas("c8", "c8", 500, 500); 
-  hClusterE->GetYaxis()->SetTitleOffset(1.4);
-  hClusterE->SetLineWidth(2);
-  hClusterE->SetLineColor(kBlue);
-  hClusterE->DrawClone();
-  c8->SaveAs("results/emcal_electrons_clusterE1.png");
-  c8->SaveAs("results/emcal_electrons_clusterE1.pdf");
-
-  TCanvas *c9 = new TCanvas("c9", "c9", 500, 500);
-  hEres->GetYaxis()->SetTitleOffset(1.4);
-  hEres->SetLineWidth(2);
-  hEres->SetLineColor(kBlue);
-  hEres->Fit("gaus");
-  hEres->DrawClone();
-  c9->SaveAs("results/emcal_electrons_Eres.png");
-  c9->SaveAs("results/emcal_electrons_Eres.pdf");
-
-  TCanvas *c10 = new TCanvas("c10", "c10", 500, 500);
-  hPthr_accepted->GetYaxis()->SetTitleOffset(1.4);
-  hPthr_accepted->SetLineWidth(2);
-  hPthr_accepted->SetLineColor(kBlue);
-  hPthr_accepted->DrawClone();
-  c10->SaveAs("results/emcal_electrons_Pthr_accepted.png");
-  c10->SaveAs("results/emcal_electrons_Pthr_accepted.pdf");
-
-  TCanvas *c11 = new TCanvas("c11", "c11", 500, 500);
-  hPratio->GetYaxis()->SetTitleOffset(1.4);
-  hPratio->SetLineWidth(2);
-  hPratio->SetLineColor(kBlue);
-  hPratio->DrawClone();
-  c11->SaveAs("results/emcal_electrons_Pratio.png");
-  c11->SaveAs("results/emcal_electrons_Pratio.pdf");
-
-  TCanvas *c12 = new TCanvas("c12", "c12", 500, 500);
-  TH1D *hPacceptance = new TH1D("hPacceptance","Ratio p_{rec}/p_{thr}", 100, -0.5, 30.5);
-  hPacceptance = (TH1D*) hPthr_accepted->Clone();
-  hPacceptance->Divide(hPthr.GetPtr());
-  hPacceptance->SetTitle("Ratio p_{rec}/p_{thr}");
-  hPacceptance->SetLineColor(kBlue);
-  hPacceptance->SetLineWidth(2);
-  hPacceptance->GetXaxis()->SetTitle("p [GeV]");
-  hPacceptance->GetYaxis()->SetTitle("p_{rec}/p_{thr}");
-  hPacceptance->GetYaxis()->SetTitleOffset(1.4);
-  hPacceptance->DrawClone();
-  c12->SaveAs("results/emcal_electrons_Pacceptance.png");
-  c12->SaveAs("results/emcal_electrons_Pacceptance.pdf");
-
-  TCanvas *c13 = new TCanvas("c13", "c13", 500, 500);
-  hEres_cut->GetYaxis()->SetTitleOffset(1.4);
-  hEres_cut->SetLineWidth(2);
-  hEres_cut->SetLineColor(kBlue);
-  hEres_cut->Fit("gaus");
-  hEres_cut->DrawClone();
-  c13->SaveAs("results/emcal_electrons_Eres_cut.png");
-  c13->SaveAs("results/emcal_electrons_Eres_cut.pdf");
-
-  TCanvas *c14 = new TCanvas("c14", "c14", 500, 500);
-  hTheta_cut->GetYaxis()->SetTitleOffset(1.4);
-  hTheta_cut->SetLineWidth(2);
-  hTheta_cut->SetLineColor(kBlue);
-  hTheta_cut->DrawClone();
-  c14->SaveAs("results/emcal_electrons_theta_cut.png"); 
-  c14->SaveAs("results/emcal_electrons_theta_cut.pdf");
-
-  TCanvas *c15 = new TCanvas("c15", "c15", 500, 500);
-  hEta_cut->GetYaxis()->SetTitleOffset(1.4);
-  hEta_cut->SetLineWidth(2);
-  hEta_cut->SetLineColor(kBlue);
-  hEta_cut->DrawClone();
-  c15->SaveAs("results/emcal_electrons_eta_cut.png");
-  c15->SaveAs("results/emcal_electrons_eta_cut.pdf");
-
-  TCanvas *c16 = new TCanvas("c16", "c16", 500, 500);
-  hEres_lmom->GetYaxis()->SetTitleOffset(1.4);
-  hEres_lmom->SetLineWidth(2);
-  hEres_lmom->SetLineColor(kBlue);
-  hEres_lmom->Fit("gaus");
-  hEres_lmom->DrawClone();
-  c16->SaveAs("results/emcal_electrons_Eres_lmom.png");
-  c16->SaveAs("results/emcal_electrons_Eres_lmom.pdf");
-
-  TCanvas *c17 = new TCanvas("c17", "c17", 500, 500);
-  hEres_hmom->GetYaxis()->SetTitleOffset(1.4);
-  hEres_hmom->SetLineWidth(2);
-  hEres_hmom->SetLineColor(kBlue);
-  hEres_hmom->Fit("gaus");
-  hEres_hmom->DrawClone();
-  c17->SaveAs("results/emcal_electrons_Eres_hmom.png");
-  c17->SaveAs("results/emcal_electrons_Eres_hmom.pdf");
-
-  TCanvas *c18 = new TCanvas("c18", "c18", 500, 500);
-  hEres_cut_lmom->GetYaxis()->SetTitleOffset(1.4);
-  hEres_cut_lmom->SetLineWidth(2);
-  hEres_cut_lmom->SetLineColor(kBlue);
-  hEres_cut_lmom->Fit("gaus");
-  hEres_cut_lmom->DrawClone();
-  c18->SaveAs("results/emcal_electrons_Eres_cut_lmom.png");
-  c18->SaveAs("results/emcal_electrons_Eres_cut_lmom.pdf");
-
-  TCanvas *c19 = new TCanvas("c19", "c19", 500, 500);
-  hEres_cut_hmom->GetYaxis()->SetTitleOffset(1.4);
-  hEres_cut_hmom->SetLineWidth(2);
-  hEres_cut_hmom->SetLineColor(kBlue);
-  hEres_cut_hmom->Fit("gaus");
-  hEres_cut_hmom->DrawClone();
-  c19->SaveAs("results/emcal_electrons_Eres_cut_hmom.png");
-  c19->SaveAs("results/emcal_electrons_Eres_cut_hmom.pdf");
-
-  TCanvas *c20 = new TCanvas("c20", "c20", 500, 500);
-  hPrec_Pthr->GetYaxis()->SetTitleOffset(1.4);
-  hPrec_Pthr->SetLineWidth(2);
-  hPrec_Pthr->SetLineColor(kBlue);
-  hPrec_Pthr->DrawClone("colz");
-  c20->SaveAs("results/emcal_electrons_Prec_Pthr.png");
-  c20->SaveAs("results/emcal_electrons_Prec_Pthr.pdf");
-
-  TCanvas *c21 = new TCanvas("c21", "c21", 500, 500);
-  hEres_Ethr->GetYaxis()->SetTitleOffset(1.4);
-  hEres_Ethr->SetLineWidth(2);
-  hEres_Ethr->SetLineColor(kBlue);
-  hEres_Ethr->DrawClone("colz");
-  c21->SaveAs("results/emcal_electrons_Eres_Ethr.png");
-  c21->SaveAs("results/emcal_electrons_Eres_Ethr.pdf");
-
-  TCanvas *c22 = new TCanvas("c22", "c22", 500, 500);
-  hEres_leta->GetYaxis()->SetTitleOffset(1.4);
-  hEres_leta->SetLineWidth(2);
-  hEres_leta->SetLineColor(kBlue);
-  hEres_leta->Fit("gaus");
-  hEres_leta->DrawClone();
-  c22->SaveAs("results/emcal_electrons_Eres_leta.png");
-  c22->SaveAs("results/emcal_electrons_Eres_leta.pdf");
-
-  TCanvas *c23 = new TCanvas("c23", "c23", 500, 500);
-  hEres_heta->GetYaxis()->SetTitleOffset(1.4);
-  hEres_heta->SetLineWidth(2);
-  hEres_heta->SetLineColor(kBlue);
-  hEres_heta->Fit("gaus");
-  hEres_heta->DrawClone();
-  c23->SaveAs("results/emcal_electrons_Eres_heta.png");
-  c23->SaveAs("results/emcal_electrons_Eres_heta.pdf");
-
-  TCanvas *c24 = new TCanvas("c24", "c24", 500, 500);
-  hEres_cut_leta->GetYaxis()->SetTitleOffset(1.4);
-  hEres_cut_leta->SetLineWidth(2);
-  hEres_cut_leta->SetLineColor(kBlue);
-  hEres_cut_leta->Fit("gaus");
-  hEres_cut_leta->DrawClone();
-  c24->SaveAs("results/emcal_electrons_Eres_cut_leta.png");
-  c24->SaveAs("results/emcal_electrons_Eres_cut_leta.pdf");
-
-  TCanvas *c25 = new TCanvas("c25", "c25", 500, 500);
-  hEres_cut_heta->GetYaxis()->SetTitleOffset(1.4);
-  hEres_cut_heta->SetLineWidth(2);
-  hEres_cut_heta->SetLineColor(kBlue);
-  hEres_cut_heta->Fit("gaus");
-  hEres_cut_heta->DrawClone();
-  c25->SaveAs("results/emcal_electrons_Eres_cut_heta.png");
-  c25->SaveAs("results/emcal_electrons_Eres_cut_heta.pdf");
-
-  TCanvas *c26 = new TCanvas("c26", "c26", 500, 500);
-  hPrec_Pthr_cut->GetYaxis()->SetTitleOffset(1.4);
-  hPrec_Pthr_cut->SetLineWidth(2);
-  hPrec_Pthr_cut->SetLineColor(kBlue);
-  hPrec_Pthr_cut->DrawClone("colz");
-  c26->SaveAs("results/emcal_electrons_Prec_Pthr_cut.png");
-  c26->SaveAs("results/emcal_electrons_Prec_Pthr_cut.pdf");
-
-  TCanvas *c27 = new TCanvas("c27", "c27", 500, 500);
-  hEres_Ethr_cut->GetYaxis()->SetTitleOffset(1.4);
-  hEres_Ethr_cut->SetLineWidth(2);
-  hEres_Ethr_cut->SetLineColor(kBlue);
-  hEres_Ethr_cut->DrawClone("colz");
-  c27->SaveAs("results/emcal_electrons_Eres_Ethr_cut.png");
-  c27->SaveAs("results/emcal_electrons_Eres_Ethr_cut.pdf");
-
-}
diff --git a/benchmarks/ecal/scripts/deprecated/emcal_electrons_reader.cxx b/benchmarks/ecal/scripts/deprecated/emcal_electrons_reader.cxx
deleted file mode 100644
index 5b47d016..00000000
--- a/benchmarks/ecal/scripts/deprecated/emcal_electrons_reader.cxx
+++ /dev/null
@@ -1,126 +0,0 @@
-//////////////////////////
-// Crystal EMCAL detector 
-// Electron dataset
-// J.KIM 07/20/2020
-//////////////////////////
-#include "HepMC3/GenEvent.h"
-#include "HepMC3/ReaderAscii.h"
-#include "HepMC3/WriterAscii.h"
-#include "HepMC3/Print.h"
-
-#include "TH1F.h"
-#include <iostream>
-#include "TStyle.h"
-
-using namespace HepMC3;
-
-void emcal_electrons_reader(double e_start = 0.0, double e_end = 30.0, const char* in_fname = "./data/emcal_electron_1000kEvt.hepmc"){
-
-      	// 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.17);
-
-	ReaderAscii hepmc_input(in_fname);
-	int         events_parsed = 0;
-	GenEvent    evt(Units::GEV, Units::MM);
-
-	// Histograms
-	TH1F* h_electrons_energy = new TH1F("h_electron_energy",  "electron energy;E [GeV];Events",         100,-0.5,30.5);
-	TH1F* h_electrons_eta    = new TH1F("h_electron_eta",     "electron #eta;#eta;Events",              100,-10.0,10.0);
-	TH1F* h_electrons_theta  = new TH1F("h_electron_theta",   "electron #theta;#theta [degree];Events", 100,-0.5,180.5);
-	TH1F* h_electrons_phi    = new TH1F("h_electron_phi",     "electron #phi;#phi [degree];Events",     100,-180.5,180.5);
-	TH2F* h_electrons_pzpt   = new TH2F("h_electrons_pzpt",   "electron pt vs pz;pt [GeV];pz [GeV]",    100,-0.5,30.5,100,-30.5,30.5);
-	TH2F* h_electrons_pxpy   = new TH2F("h_electrons_pxpy",   "electron px vs py;px [GeV];py [GeV]",    100,-30.5,30.5,100,-30.5,30.5);
-	TH3F* h_electrons_p      = new TH3F("h_electron_p",       "electron p;px [GeV];py [GeV];pz [GeV]",  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;
-
-    		for(const auto& v : evt.vertices() ) {
-      			for(const auto& p : v->particles_out() ) {
-        			if(p->pid() == 11) {
-          			h_electrons_energy->Fill(p->momentum().e());
-				h_electrons_eta->Fill(p->momentum().eta());
-				h_electrons_theta->Fill(p->momentum().theta()*TMath::RadToDeg());
-				h_electrons_phi->Fill(p->momentum().phi()*TMath::RadToDeg());
-				h_electrons_pzpt->Fill(TMath::Sqrt(p->momentum().px()*p->momentum().px() + p->momentum().py()*p->momentum().py()),p->momentum().pz());
-				h_electrons_pxpy->Fill(p->momentum().px(),p->momentum().py());
-				h_electrons_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_electrons_energy->GetYaxis()->SetTitleOffset(1.8);
-	h_electrons_energy->SetLineWidth(2);
-	h_electrons_energy->SetLineColor(kBlue);
-  	h_electrons_energy->DrawClone();
-  	c->SaveAs("results/input_emcal_electrons_energy.png");
-  	c->SaveAs("results/input_emcal_electrons_energy.pdf");
-
-        TCanvas* c1 = new TCanvas("c1","c1", 500, 500);
-	h_electrons_eta->GetYaxis()->SetTitleOffset(1.9);
-	h_electrons_eta->SetLineWidth(2);
-	h_electrons_eta->SetLineColor(kBlue);
-        h_electrons_eta->DrawClone();
-        c1->SaveAs("results/input_emcal_electrons_eta.png");
-        c1->SaveAs("results/input_emcal_electrons_eta.pdf");
-
-        TCanvas* c2 = new TCanvas("c2","c2", 500, 500);
-	h_electrons_theta->GetYaxis()->SetTitleOffset(1.8);
-	h_electrons_theta->SetLineWidth(2);
-	h_electrons_theta->SetLineColor(kBlue);
-        h_electrons_theta->DrawClone();
-	c2->SaveAs("results/input_emcal_electrons_theta.png");
-	c2->SaveAs("results/input_emcal_electrons_theta.pdf");
-
-        TCanvas* c3 = new TCanvas("c3","c3", 500, 500);
-	h_electrons_phi->GetYaxis()->SetTitleOffset(1.8);
-	h_electrons_phi->SetLineWidth(2);
-	h_electrons_phi->GetYaxis()->SetRangeUser(0.0,h_electrons_phi->GetMaximum()+300.0);
-	h_electrons_phi->SetLineColor(kBlue);
-        h_electrons_phi->DrawClone();
-	c3->SaveAs("results/input_emcal_electrons_phi.png");
-	c3->SaveAs("results/input_emcal_electrons_phi.pdf");
-
-	TCanvas* c4 = new TCanvas("c4","c4", 500, 500);
-	h_electrons_pzpt->GetYaxis()->SetTitleOffset(1.4);
-	h_electrons_pzpt->SetLineWidth(2);
-	h_electrons_pzpt->SetLineColor(kBlue);
-	h_electrons_pzpt->DrawClone("COLZ");
-	c4->SaveAs("results/input_emcal_electrons_pzpt.png");
-	c4->SaveAs("results/input_emcal_electrons_pzpt.pdf");
-
-	TCanvas* c5 = new TCanvas("c5","c5", 500, 500);
-	h_electrons_pxpy->GetYaxis()->SetTitleOffset(1.4);
-	h_electrons_pxpy->SetLineWidth(2);
-	h_electrons_pxpy->SetLineColor(kBlue);
-	h_electrons_pxpy->DrawClone("COLZ");
-	c5->SaveAs("results/input_emcal_electrons_pxpy.png");
-	c5->SaveAs("results/input_emcal_electrons_pxpy.pdf");
-
- 	TCanvas* c6 = new TCanvas("c6","c6", 500, 500);
-	h_electrons_p->GetYaxis()->SetTitleOffset(1.8);
-	h_electrons_p->GetXaxis()->SetTitleOffset(1.6);
-	h_electrons_p->GetZaxis()->SetTitleOffset(1.6);
-	h_electrons_p->SetLineWidth(2);
-	h_electrons_p->SetLineColor(kBlue);
-        h_electrons_p->DrawClone();
-	c6->SaveAs("results/input_emcal_electrons_p.png");
-	c6->SaveAs("results/input_emcal_electrons_p.pdf");
-
-}
diff --git a/benchmarks/ecal/scripts/deprecated/emcal_pi0.cxx b/benchmarks/ecal/scripts/deprecated/emcal_pi0.cxx
deleted file mode 100644
index 3d6fd0c5..00000000
--- a/benchmarks/ecal/scripts/deprecated/emcal_pi0.cxx
+++ /dev/null
@@ -1,87 +0,0 @@
-//////////////////////////////////////////////////////////////
-// Crystal EMCAL detector 
-// Single Pion dataset
-// J.KIM 10/18/2020
-//////////////////////////////////////////////////////////////
-#include "HepMC3/GenEvent.h"
-#include "HepMC3/ReaderAscii.h"
-#include "HepMC3/WriterAscii.h"
-#include "HepMC3/Print.h"
-
-#include <iostream>
-#include<random>
-#include<cmath>
-#include <math.h>
-#include <TMath.h>
-
-using namespace HepMC3;
-
-void emcal_pi0(int n_events = 1e6, const char* out_fname = "./data/emcal_pi0_0GeVto30GeV_1000kEvt.hepmc")
-{
-  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 * (160.0 / 180.0));
-  double cos_theta_max = std::cos(M_PI * (175.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 111 - pi0
-    // 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
-    // Momentum starting at 0 GeV/c to 30 GeV/c
-    Double_t p        = r1->Uniform(0.0, 30.0);
-    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);
-
-    // Generates random vectors, uniformly distributed over the surface of a
-    // sphere of given radius, in this case momentum.
-    //r1->Sphere(px, py, pz, p);
-
-    // type 1 is final state
-    // pdgid 111 - pi0 135 MeV/c^2
-    GenParticlePtr p3 = std::make_shared<GenParticle>(
-        FourVector(
-            px, py, pz,
-            sqrt(p*p + (0.134976*0.134976))),
-        111, 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/ecal/scripts/deprecated/emcal_pi0_reader.cxx b/benchmarks/ecal/scripts/deprecated/emcal_pi0_reader.cxx
deleted file mode 100644
index 95f631df..00000000
--- a/benchmarks/ecal/scripts/deprecated/emcal_pi0_reader.cxx
+++ /dev/null
@@ -1,138 +0,0 @@
-//////////////////////////
-// Crystal EMCAL detector 
-// Singe Pion dataset
-// J.KIM 10/18/2020
-//////////////////////////
-#include "HepMC3/GenEvent.h"
-#include "HepMC3/ReaderAscii.h"
-#include "HepMC3/WriterAscii.h"
-#include "HepMC3/Print.h"
-
-#include "TH1F.h"
-#include <iostream>
-#include "TStyle.h"
-
-using namespace HepMC3;
-
-void emcal_pi0_reader(const char* in_fname = "./data/emcal_pi0_0GeVto30GeV_1000kEvt.hepmc"){
-
-      	// 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.17);
-
-	ReaderAscii hepmc_input(in_fname);
-	int         events_parsed = 0;
-	GenEvent    evt(Units::GEV, Units::MM);
-
-	// Histograms
-	TH1F* h_pi0_energy = new TH1F("h_pi0_energy",  "pi0 energy;E [GeV];Events",         100,-0.5,30.5);
-	TH1F* h_pi0_eta    = new TH1F("h_pi0_eta",     "pi0 #eta;#eta;Events",              100,-10.0,10.0);
-	TH1F* h_pi0_theta  = new TH1F("h_pi0_theta",   "pi0 #theta;#theta [degree];Events", 100,-0.5,180.5);
-	TH1F* h_pi0_phi    = new TH1F("h_pi0_phi",     "pi0 #phi;#phi [degree];Events",     100,-180.5,180.5);
-	TH2F* h_pi0_pzpt   = new TH2F("h_pi0_pzpt",    "pi0 pt vs pz;pt [GeV];pz [GeV]",    100,-0.5,30.5,100,-30.5,30.5);
-	TH2F* h_pi0_pxpy   = new TH2F("h_pi0_pxpy",    "pi0 px vs py;px [GeV];py [GeV]",    100,-30.5,30.5,100,-30.5,30.5);
-	TH3F* h_pi0_p      = new TH3F("h_pi0_p",       "pi0 p;px [GeV];py [GeV];pz [GeV]",  100,-30.5,30.5,100,-30.5,30.5,100,-30.5,30.5);
-	TH1F* h_pi0_mom    = new TH1F("h_pi0_mom",     "pi0 momentum;p [GeV];Events",       100,-0.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;
-
-    		for(const auto& v : evt.vertices() ) {
-      			for(const auto& p : v->particles_out() ) {
-        			if(p->pid() == 111) {
-          			h_pi0_energy->Fill(p->momentum().e());
-				h_pi0_eta->Fill(p->momentum().eta());
-				h_pi0_theta->Fill(p->momentum().theta()*TMath::RadToDeg());
-				h_pi0_phi->Fill(p->momentum().phi()*TMath::RadToDeg());
-				h_pi0_pzpt->Fill(TMath::Sqrt(p->momentum().px()*p->momentum().px() + p->momentum().py()*p->momentum().py()),p->momentum().pz());
-				h_pi0_pxpy->Fill(p->momentum().px(),p->momentum().py());
-				h_pi0_p->Fill(p->momentum().px(),p->momentum().py(),p->momentum().pz());
-				h_pi0_mom->Fill(TMath::Sqrt(p->momentum().px()*p->momentum().px() + 
-							p->momentum().py()*p->momentum().py() + 
-							p->momentum().pz()*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_pi0_energy->GetYaxis()->SetTitleOffset(1.8);
-	h_pi0_energy->SetLineWidth(2);
-	h_pi0_energy->SetLineColor(kBlue);
-  	h_pi0_energy->DrawClone();
-  	c->SaveAs("results/input_emcal_pi0_energy_0GeVto30GeV.png");
-  	c->SaveAs("results/input_emcal_pi0_energy_0GeVto30GeV.pdf");
-
-        TCanvas* c1 = new TCanvas("c1","c1", 500, 500);
-	h_pi0_eta->GetYaxis()->SetTitleOffset(1.9);
-	h_pi0_eta->SetLineWidth(2);
-	h_pi0_eta->SetLineColor(kBlue);
-        h_pi0_eta->DrawClone();
-        c1->SaveAs("results/input_emcal_pi0_eta_0GeVto30GeV.png");
-        c1->SaveAs("results/input_emcal_pi0_eta_0GeVto30GeV.pdf");
-
-        TCanvas* c2 = new TCanvas("c2","c2", 500, 500);
-	h_pi0_theta->GetYaxis()->SetTitleOffset(1.8);
-	h_pi0_theta->SetLineWidth(2);
-	h_pi0_theta->SetLineColor(kBlue);
-        h_pi0_theta->DrawClone();
-	c2->SaveAs("results/input_emcal_pi0_theta_0GeVto30GeV.png");
-	c2->SaveAs("results/input_emcal_pi0_theta_0GeVto30GeV.pdf");
-
-        TCanvas* c3 = new TCanvas("c3","c3", 500, 500);
-	h_pi0_phi->GetYaxis()->SetTitleOffset(1.8);
-	h_pi0_phi->SetLineWidth(2);
-	h_pi0_phi->GetYaxis()->SetRangeUser(0.0,h_pi0_phi->GetMaximum()+300.0);
-	h_pi0_phi->SetLineColor(kBlue);
-        h_pi0_phi->DrawClone();
-	c3->SaveAs("results/input_emcal_pi0_phi_0GeVto30GeV.png");
-	c3->SaveAs("results/input_emcal_pi0_phi_0GeVto30GeV.pdf");
-
-	TCanvas* c4 = new TCanvas("c4","c4", 500, 500);
-	h_pi0_pzpt->GetYaxis()->SetTitleOffset(1.4);
-	h_pi0_pzpt->SetLineWidth(2);
-	h_pi0_pzpt->SetLineColor(kBlue);
-	h_pi0_pzpt->DrawClone("COLZ");
-	c4->SaveAs("results/input_emcal_pi0_pzpt_0GeVto30GeV.png");
-	c4->SaveAs("results/input_emcal_pi0_pzpt_0GeVto30GeV.pdf");
-
-	TCanvas* c5 = new TCanvas("c5","c5", 500, 500);
-	h_pi0_pxpy->GetYaxis()->SetTitleOffset(1.4);
-	h_pi0_pxpy->SetLineWidth(2);
-	h_pi0_pxpy->SetLineColor(kBlue);
-	h_pi0_pxpy->DrawClone("COLZ");
-	c5->SaveAs("results/input_emcal_pi0_pxpy_0GeVto30GeV.png");
-	c5->SaveAs("results/input_emcal_pi0_pxpy_0GeVto30GeV.pdf");
-
- 	TCanvas* c6 = new TCanvas("c6","c6", 500, 500);
-	h_pi0_p->GetYaxis()->SetTitleOffset(1.8);
-	h_pi0_p->GetXaxis()->SetTitleOffset(1.6);
-	h_pi0_p->GetZaxis()->SetTitleOffset(1.6);
-	h_pi0_p->SetLineWidth(2);
-	h_pi0_p->SetLineColor(kBlue);
-        h_pi0_p->DrawClone();
-	c6->SaveAs("results/input_emcal_pi0_p_0GeVto30GeV.png");
-	c6->SaveAs("results/input_emcal_pi0_p_0GeVto30GeV.pdf");
-
-        TCanvas* c7 = new TCanvas("c7","c7", 500, 500);
-        h_pi0_mom->GetYaxis()->SetTitleOffset(1.8);
-        h_pi0_mom->SetLineWidth(2);
-        h_pi0_mom->SetLineColor(kBlue);
-        h_pi0_mom->DrawClone();
-        c7->SaveAs("results/input_emcal_pi0_mom_0GeVto30GeV.png");
-        c7->SaveAs("results/input_emcal_pi0_mom_0GeVto30GeV.pdf");
-
-}
diff --git a/benchmarks/ecal/scripts/deprecated/emcal_pion_anlaysis.cxx b/benchmarks/ecal/scripts/deprecated/emcal_pion_anlaysis.cxx
deleted file mode 100644
index ca1c7fa2..00000000
--- a/benchmarks/ecal/scripts/deprecated/emcal_pion_anlaysis.cxx
+++ /dev/null
@@ -1,251 +0,0 @@
-////////////////////////////////////////
-// Read reconstruction ROOT output file
-// Plot variables
-////////////////////////////////////////
-
-#include "ROOT/RDataFrame.hxx"
-#include <iostream>
-
-#include "dd4pod/Geant4ParticleCollection.h"
-#include "dd4pod/CalorimeterHitCollection.h"
-#include "eicd/CalorimeterHitCollection.h"
-#include "eicd/CalorimeterHitData.h"
-#include "eicd/ClusterCollection.h"
-#include "eicd/ClusterData.h"
-
-#include "TCanvas.h"
-#include "TStyle.h"
-#include "TMath.h"
-#include "TH1.h"
-#include "TH1D.h"
-
-// If root cannot find a dd4pod header make sure to add the include dir to ROOT_INCLUDE_PATH:
-//  export ROOT_INCLUDE_PATH=$HOME/include:$ROOT_INCLUDE_PATH
-//  export ROOT_INCLUDE_PATH=/home/jihee/stow/development/include:$ROOT_INCLUDE_PATH
-
-using ROOT::RDataFrame;
-using namespace ROOT::VecOps;
-
-void emcal_pion_anlaysis(const char* input_fname = "rec_emcal_uniform_pi0s.root")
-{
-  // 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();
-  ROOT::RDataFrame d0("events", input_fname);
-
-  // Number of Clusters
-  auto ncluster = [] (const std::vector<eic::ClusterData>& evt) {
-	  //std::cout << (int) evt.size() << endl;
-	  return (int) evt.size(); };
-
-  // Angle between two photons [degree]
-  auto theta_two_photons = [] (const std::vector<eic::ClusterData>& evt) {
-        std::vector<double> result;
-        double cluster_posx[2];
-        double cluster_posy[2];
-        double cluster_posz[2];
-        double cluster_energy[2];
-
-        if((int) evt.size() == 2) {
-                std::cout << "Found it" << endl;
-                auto n = 0;
-                for(const auto& i: evt) {
-                        //std::cout << n << ": " <<  i.energy << endl;
-                        cluster_posx[n] = i.position.x;
-                        cluster_posy[n] = i.position.y;
-                        cluster_posz[n] = i.position.z;
-                        cluster_energy[n] = i.energy;
-                        n++;
-                }
-        }
-        // Calculate invariant mass
-        auto dot_product_pos_clusters = cluster_posx[0]*cluster_posx[1] + cluster_posy[0]*cluster_posy[1] + cluster_posz[0]*cluster_posz[1];
-        auto mag_pos2_cluster_1       = (cluster_posx[0]*cluster_posx[0]) + (cluster_posy[0]*cluster_posy[0]) + (cluster_posz[0]*cluster_posz[0]);
-        auto mag_pos2_cluster_2       = (cluster_posx[1]*cluster_posx[1]) + (cluster_posy[1]*cluster_posy[1]) + (cluster_posz[1]*cluster_posz[1]);
-        auto cosine_clusters          = (dot_product_pos_clusters/TMath::Sqrt(mag_pos2_cluster_1*mag_pos2_cluster_2));
-        auto theta_photons            = TMath::ACos(cosine_clusters)*TMath::RadToDeg();
-        result.push_back(theta_photons);
-  return result;
-  };
-
-  // Invariant Mass [MeV]
-  auto inv_mass = [] (const std::vector<eic::ClusterData>& evt) {
-  	std::vector<double> result;
-	double cluster_posx[2];
-        double cluster_posy[2];
-        double cluster_posz[2];
-        double cluster_energy[2];
-
-	if((int) evt.size() == 2) {
-		//std::cout << "Found it" << endl;
-		auto n = 0;
-        	for(const auto& i: evt) {
-			//std::cout << n << ": " <<  i.energy << endl;
-			cluster_posx[n] = i.position.x;
-			cluster_posy[n] = i.position.y;
-			cluster_posz[n] = i.position.z;
-			cluster_energy[n] = i.energy;
-			n++;
-        	}
-	}
-	// Calculate invariant mass
-	auto dot_product_pos_clusters = cluster_posx[0]*cluster_posx[1] + cluster_posy[0]*cluster_posy[1] + cluster_posz[0]*cluster_posz[1];
-        auto mag_pos2_cluster_1       = (cluster_posx[0]*cluster_posx[0]) + (cluster_posy[0]*cluster_posy[0]) + (cluster_posz[0]*cluster_posz[0]);
-        auto mag_pos2_cluster_2       = (cluster_posx[1]*cluster_posx[1]) + (cluster_posy[1]*cluster_posy[1]) + (cluster_posz[1]*cluster_posz[1]);
-        auto cosine_clusters          = (dot_product_pos_clusters/TMath::Sqrt(mag_pos2_cluster_1*mag_pos2_cluster_2));
-        auto theta_photons            = TMath::ACos(cosine_clusters)*TMath::RadToDeg();
-        auto invariant_mass           = TMath::Sqrt(2.0*cluster_energy[0]*cluster_energy[1]*(1.0 - cosine_clusters))*1.e+3;
-	result.push_back(invariant_mass);
-  return result;
-  };
-
-  // Thrown Energy [GeV]
-  auto E_thr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
-        std::vector<double> result;
-        result.push_back(TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass));
-  return result;
-  };
-
-  // Cluster Energy [GeV]
-  auto clusterE = [] (const std::vector<eic::ClusterData>& evt) {
-        std::vector<double> result;
-        auto total_eng = 0.0;
-	for (const auto& i: evt) {
-        	 total_eng += i.energy;
-	}
-	result.push_back(total_eng);
-  return result;
-  };
-
-  // Energy Resolution
-  auto E_res = [] (const std::vector<double>& Erec, const std::vector<double>& Ethr) {
-        std::vector<double> result;
-        for (const auto& E1 : Ethr) {
-                for (const auto& E2 : Erec) {
-                        result.push_back((E2-E1)/E1);
-                }
-        }
-  return result;
-  };  
-
-  // Define variables
-  auto d1 = d0.Define("ncluster", ncluster, {"EcalClusters"})
-	      .Define("theta_two_photons", theta_two_photons, {"EcalClusters"})
-	      .Define("inv_mass", inv_mass, {"EcalClusters"})
-	      .Define("E_thr",    E_thr,    {"mcparticles2"})
-	      .Define("clusterE", clusterE, {"EcalClusters"})
-	      .Define("E_res",    E_res,    {"clusterE","E_thr"})
-	      ;
-
-  // Define Histograms
-  auto hNCluster = d1.Histo1D({"hNCluster", "Number of Clusters; # of Clusters; Events",    5, -0.5, 4.5},      "ncluster");
-  auto hInvMass  = d1.Histo1D({"hInvMass",  "Invariant Mass; mass [MeV]; Events", 60,0.0,300.0}, "inv_mass");
-
-  // Select Events with One Cluster
-  auto d2 = d1.Filter("ncluster==2");
-  auto hInvMass_NC2  = d2.Histo1D({"hInvMass_NC2",  "Invariant Mass; mass [MeV]; Events", 60,0.0,300.0}, "inv_mass");
-  auto hThetaTwoPhotons_NC2  = d2.Histo1D({"hThetaTwoPhotons_NC2",  "Angle between Two Photons; angle [degree]; Events", 100,0.0,30.0}, "theta_two_photons");
-  auto hEres_NC2         = d2.Histo1D({"hEres_NC2",      "Energy Resolution; #DeltaE/E; Events",             100, -0.5, 0.5},  "E_res"); 
-
-  // Select Events with Edge CUT
-
-  auto d3 = d2.Filter([=] (const std::vector<eic::ClusterData>& evt) {
-        double cluster_posx[2];
-        double cluster_posy[2];
-
-        auto n = 0;
-        for(const auto& i: evt) {
-        	cluster_posx[n] = i.position.x;
-                cluster_posy[n] = i.position.y;
-                n++;
-        }
-	auto radial_dist_c1 = TMath::Sqrt(cluster_posx[0]*cluster_posx[0] + cluster_posy[0]*cluster_posy[0]);
-	auto radial_dist_c2 = TMath::Sqrt(cluster_posx[1]*cluster_posx[1] + cluster_posy[1]*cluster_posy[1]);
-
-        if (radial_dist_c1 > 18.0 && radial_dist_c1 < 30.0 && radial_dist_c2 > 18.0 && radial_dist_c2 < 30.0)
-		return true;
-        return false;}, {"EcalClusters"});
-  auto hEres_NC2_CUT         = d3.Histo1D({"hEres_NC2_CUT",      "Energy Resolution; #DeltaE/E; Events",             100, -0.5, 0.5},  "E_res"); 
-  auto hInvMass_NC2_CUT  = d3.Histo1D({"hInvMass_NC2_CUT",  "Invariant Mass; mass [MeV]; Events", 60,0.0,300.0}, "inv_mass");
-
-  // Event Counts
-  auto nevents_thrown    = d1.Count();
-  auto nevents_cluster1  = d2.Count();
-  auto nevents_cluster1cut = d3.Count();
-
-  std::cout << "Number of Events: " << (*nevents_cluster1) << " / " << (*nevents_thrown) << " = single cluster events/all \n";
-  std::cout << "Number of Events: " << (*nevents_cluster1cut) << " / " << (*nevents_thrown) << " = single cluster events that passed edge cut/all \n";
-
-  // Draw Histograms
-  TCanvas *c1 = new TCanvas("c1", "c1", 500, 500);
-  c1->SetLogy(1); 
-  hNCluster->GetYaxis()->SetTitleOffset(1.6);
-  hNCluster->SetLineWidth(2);
-  hNCluster->SetLineColor(kBlue);
-  hNCluster->DrawClone();
-  c1->SaveAs("results/emcal_pi0s_ncluster.png");
-  c1->SaveAs("results/emcal_pi0s_ncluster.pdf");
-
-  TCanvas *c2 = new TCanvas("c2", "c2", 500, 500); 
-  hInvMass->GetYaxis()->SetTitleOffset(1.4);
-  hInvMass->SetLineWidth(2);
-  hInvMass->SetLineColor(kBlue);
-  hInvMass->Fit("gaus");
-  hInvMass->DrawClone();
-  c2->SaveAs("results/emcal_pi0s_invariant_mass.png");
-  c2->SaveAs("results/emcal_pi0s_invariant_mass.pdf");
-
-  TCanvas *c3 = new TCanvas("c3", "c3", 500, 500); 
-  hInvMass_NC2->GetYaxis()->SetTitleOffset(1.4);
-  hInvMass_NC2->SetLineWidth(2);
-  hInvMass_NC2->SetLineColor(kBlue);
-  hInvMass_NC2->Fit("gaus");
-  hInvMass_NC2->DrawClone();
-  c3->SaveAs("results/emcal_pi0s_invariant_mass_nc2.png");
-  c3->SaveAs("results/emcal_pi0s_invariant_mass_nc2.pdf");
-
-  TCanvas *c4 = new TCanvas("c4", "c4", 500, 500);
-  hThetaTwoPhotons_NC2->GetYaxis()->SetTitleOffset(1.4);
-  hThetaTwoPhotons_NC2->SetLineWidth(2);
-  hThetaTwoPhotons_NC2->SetLineColor(kBlue);
-  hThetaTwoPhotons_NC2->DrawClone();
-  c4->SaveAs("results/emcal_pi0s_angle_two_photons_nc2.png");
-  c4->SaveAs("results/emcal_pi0s_angle_two_photons_nc2.pdf");
-
-  TCanvas *c5 = new TCanvas("c5", "c5", 500, 500);
-  hEres_NC2->GetYaxis()->SetTitleOffset(1.4);
-  hEres_NC2->SetLineWidth(2);
-  hEres_NC2->SetLineColor(kBlue);
-  hEres_NC2->Fit("gaus");
-  hEres_NC2->DrawClone();
-  c5->SaveAs("results/emcal_pi0s_Eres_nc2.png");
-  c5->SaveAs("results/emcal_pi0s_Eres_nc2.pdf");
- 
-  TCanvas *c6 = new TCanvas("c6", "c6", 500, 500);
-  hEres_NC2_CUT->GetYaxis()->SetTitleOffset(1.4);
-  hEres_NC2_CUT->SetLineWidth(2);
-  hEres_NC2_CUT->SetLineColor(kBlue);
-  hEres_NC2_CUT->Fit("gaus");
-  hEres_NC2_CUT->DrawClone();
-  c6->SaveAs("results/emcal_pi0s_Eres_nc2_cut.png");
-  c6->SaveAs("results/emcal_pi0s_Eres_nc2_cut.pdf");
- 
-  TCanvas *c7 = new TCanvas("c7", "c7", 500, 500); 
-  hInvMass_NC2_CUT->GetYaxis()->SetTitleOffset(1.4);
-  hInvMass_NC2_CUT->SetLineWidth(2);
-  hInvMass_NC2_CUT->SetLineColor(kBlue);
-  hInvMass_NC2_CUT->Fit("gaus");
-  hInvMass_NC2_CUT->DrawClone();
-  c7->SaveAs("results/emcal_pi0s_invariant_mass_nc2_cut.png");
-  c7->SaveAs("results/emcal_pi0s_invariant_mass_nc2_cut.pdf");
- 
-}
diff --git a/benchmarks/ecal/scripts/deprecated/full_emcal_electrons.cxx b/benchmarks/ecal/scripts/deprecated/full_emcal_electrons.cxx
deleted file mode 100644
index f849c076..00000000
--- a/benchmarks/ecal/scripts/deprecated/full_emcal_electrons.cxx
+++ /dev/null
@@ -1,89 +0,0 @@
-//////////////////////////////////////////////////////////////
-// Tungsten Sampling EMCAL detector
-// Barrel and Endcap 
-// Single Electron dataset
-// J.KIM 02/12/2021
-//
-//////////////////////////////////////////////////////////////
-#include "HepMC3/GenEvent.h"
-#include "HepMC3/ReaderAscii.h"
-#include "HepMC3/WriterAscii.h"
-#include "HepMC3/Print.h"
-
-#include <iostream>
-#include<random>
-#include<cmath>
-#include <math.h>
-#include <TMath.h>
-
-using namespace HepMC3;
-
-void full_emcal_electrons(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/emcal_electron_0GeVto30GeV_1000kEvt.hepmc")
-{
-  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 111 - pi0
-    // 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
-    // Momentum starting at 0 GeV/c to 30 GeV/c
-    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);
-
-    // Generates random vectors, uniformly distributed over the surface of a
-    // sphere of given radius, in this case momentum.
-    //r1->Sphere(px, py, pz, p);
-
-    // type 1 is final state
-    // pdgid 11 - electron 0.510 MeV/c^2
-    GenParticlePtr p3 = std::make_shared<GenParticle>(
-        FourVector(
-            px, py, pz,
-            sqrt(p*p + (0.000511 * 0.000511))),
-        11, 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/ecal/scripts/deprecated/full_emcal_electrons_analysis.cxx b/benchmarks/ecal/scripts/deprecated/full_emcal_electrons_analysis.cxx
deleted file mode 100644
index 25fbecb2..00000000
--- a/benchmarks/ecal/scripts/deprecated/full_emcal_electrons_analysis.cxx
+++ /dev/null
@@ -1,70 +0,0 @@
-////////////////////////////////////////
-// Read reconstruction ROOT output file
-// Plot variables
-////////////////////////////////////////
-
-#include "ROOT/RDataFrame.hxx"
-#include <iostream>
-
-#include "dd4pod/Geant4ParticleCollection.h"
-#include "eicd/ClusterCollection.h"
-#include "eicd/ClusterData.h"
-
-#include "TCanvas.h"
-#include "TStyle.h"
-#include "TMath.h"
-#include "TH1.h"
-#include "TH1D.h"
-
-// If root cannot find a dd4pod header make sure to add the include dir to ROOT_INCLUDE_PATH:
-//  export ROOT_INCLUDE_PATH=$HOME/include:$ROOT_INCLUDE_PATH
-//  export ROOT_INCLUDE_PATH=/home/jihee/stow/development/include:$ROOT_INCLUDE_PATH
-
-using ROOT::RDataFrame;
-using namespace ROOT::VecOps;
-
-void full_emcal_electrons_analysis(const char* input_fname = "rec_emcal_uniform_electrons.root")
-{
-  // 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();
-  ROOT::RDataFrame d0("events", input_fname);
-
-  // Thrown Energy [GeV]
-  auto E_thr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
-  	std::vector<float> result;
-    	result.push_back(TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass));
-  return result;
-  };
-
-  // Define variables
-  auto d1 = d0.Define("E_thr",    E_thr,    {"mcparticles2"})
-	      ;
-
-  // Define Histograms
-  auto hEthr     = d1.Histo1D({"hEthr",     "Thrown Energy; E_{thr}[GeV]; Events",          100, -0.5, 30.5},   "E_thr");
-
-  // Event Counts
-  auto nevents_thrown    = d1.Count();
-
-  // Print out number of events
-  std::cout << "Number of Events: " << (*nevents_thrown) << " = all thrown events \n";
-
-  // Draw Histogram
-  TCanvas *c1 = new TCanvas("c1", "c1", 500, 500);
-  hEthr->GetYaxis()->SetTitleOffset(1.4);
-  hEthr->SetLineWidth(2);
-  hEthr->SetLineColor(kBlue);
-  hEthr->DrawClone();
-  c1->SaveAs("results/emcal_electrons_Ethr.png");
-  c1->SaveAs("results/emcal_electrons_Ethr.pdf");
-}
diff --git a/benchmarks/ecal/scripts/deprecated/makeplot.C b/benchmarks/ecal/scripts/deprecated/makeplot.C
deleted file mode 100644
index e7b88e26..00000000
--- a/benchmarks/ecal/scripts/deprecated/makeplot.C
+++ /dev/null
@@ -1,249 +0,0 @@
-////////////////////////////////////////
-// Read reconstruction ROOT output file
-// Plot variables
-////////////////////////////////////////
-int makeplot(double e_start = 1.0, double e_end = 1.0, const char* input_fname = "sim_output/sim_emcal_electrons_output.root", const char *fout = "results/rec_1000k_info_1GeV.txt")
-{
-  // Setting figures
-  gROOT->SetStyle("Plain");
-  gStyle->SetLineWidth(3);
-  gStyle->SetOptStat("nem");
-  gStyle->SetPadTickX(1);
-  gStyle->SetPadTickY(1);
-  gStyle->SetPadGridX(1);
-  gStyle->SetPadGridY(1);
-  gStyle->SetPadLeftMargin(0.14);
-  gStyle->SetPadRightMargin(0.14);
-
-  // Output ASCII file
-  FILE *fp = fopen(fout,"a");
-  if(!fp)
-        fprintf(stderr,"error: can't open %s\n", fout);
-
-  // Input ROOT file
-  TFile *f = new TFile(input_fname,"READ");
-  TTree *t = (TTree *)f->Get("events");
-
-  // Set Branch status and addressed
-  t->SetMakeClass(1);
-  t->SetBranchStatus("*", 0);
-
-  Int_t EcalClusters_;
-  t->SetBranchStatus("EcalClusters", 1);
-  t->SetBranchAddress("EcalClusters", &EcalClusters_);
-
-  Int_t RecoEcalHits_;
-  t->SetBranchStatus("RecoEcalHits", 1);
-  t->SetBranchAddress("RecoEcalHits", &RecoEcalHits_);
-
-  const Int_t kMaxEcalClusters = 100000;
-  Double_t cluster_x_pos[kMaxEcalClusters];
-  Double_t cluster_y_pos[kMaxEcalClusters];
-  Double_t cluster_z_pos[kMaxEcalClusters];
-  Float_t cluster_energy[kMaxEcalClusters];
-  t->SetBranchStatus("EcalClusters.position.x",1);
-  t->SetBranchStatus("EcalClusters.position.y",1);
-  t->SetBranchStatus("EcalClusters.position.z",1);
-  t->SetBranchStatus("EcalClusters.energy",1);
-  t->SetBranchAddress("EcalClusters.position.x",cluster_x_pos);
-  t->SetBranchAddress("EcalClusters.position.y",cluster_y_pos);
-  t->SetBranchAddress("EcalClusters.position.z",cluster_z_pos);
-  t->SetBranchAddress("EcalClusters.energy",cluster_energy);
-
-  const Int_t kMaxRecoEcalHits = 100000;
-  Double_t rec_x_pos[kMaxRecoEcalHits];
-  Double_t rec_y_pos[kMaxRecoEcalHits];
-  Double_t rec_energy[kMaxRecoEcalHits];
-  t->SetBranchStatus("RecoEcalHits.position.x",1);
-  t->SetBranchStatus("RecoEcalHits.position.y",1);
-  t->SetBranchStatus("RecoEcalHits.energy",1);
-  t->SetBranchAddress("RecoEcalHits.position.x",rec_x_pos);
-  t->SetBranchAddress("RecoEcalHits.position.y",rec_y_pos);
-  t->SetBranchAddress("RecoEcalHits.energy",rec_energy);
-
-  // Setting for Canvas
-  TCanvas *c1  = new TCanvas("c1", "c1",  500, 500);
-  TCanvas *c2  = new TCanvas("c2", "c2",  500, 500);
-  TCanvas *c3  = new TCanvas("c3", "c3",  500, 500);
-  TCanvas *c4  = new TCanvas("c4", "c4",  500, 500);
-  TCanvas *c5  = new TCanvas("c5", "c5",  500, 500);
-  TCanvas *c6  = new TCanvas("c6", "c6",  500, 500);
-  TCanvas *c7  = new TCanvas("c7", "c7",  500, 500);
-  TCanvas *c8  = new TCanvas("c8", "c8",  500, 500);
-  TCanvas *c9  = new TCanvas("c9", "c9",  500, 500);
-  
-  // Declare histograms
-  TH1D *h1 = new TH1D("h1","Scattering Angle(#theta)",          100,130.0,180.0);
-  TH1D *h2 = new TH1D("h2","Pseudo-rapidity(#eta)",             100,-5.0,0.0);
-  TH2D *h3 = new TH2D("h3","Cluster E vs Pseudo-rapidity",      100,e_start-0.5,e_end+0.5,100,-5.0,0.0);
-  TH1D *h4 = new TH1D("h4","Reconstructed energy per event",    100,e_start-0.5,e_end+0.5);
-  TH1D *h5 = new TH1D("h5","Number of Clusters per event",      5,-0.5,4.5);
-  TH1D *h6 = new TH1D("h6","Scattering Angle(#theta) with CUT", 100,130.0,180.0);
-  TH1D *h7 = new TH1D("h7","Pseudo-rapidity(#eta) with CUT",    100,-5.0,0.0);
-  TH2D *h8 = new TH2D("h8","Cluster Hit Position",              62,-62.0,62.0,62,-62.0,62.0);
-  TH2D *h9 = new TH2D("h9","All Hit Position",                  62,-62.0,62.0,62,-62.0,62.0);
-
-  // Declare ellipse for boundary of crystal calorimeter
-  TEllipse *ell1 = new TEllipse(0.0, 0.0, 60.0, 60.0);
-  ell1->SetFillStyle(4000);
-  TEllipse *ell2 = new TEllipse(0.0, 0.0, 12.0, 12.0);
-  ell2->SetFillStyle(4000);
-
-  // Total number of entries
-  Int_t nentries = t->GetEntries();
-
-  // Variables are used in calculation
-  Double_t r;              // Radius [cm]
-  Double_t phi;            // Azimuth [degree]
-  Double_t theta;          // Inclination [degree]
-  Double_t eta;            // Pseudo-rapidity [unitless]
-  Float_t  cluster_e;      // Cluster energy [GeV]
-  Float_t  total_cluster_e;// Add up clusters per event [GeV]
-  Double_t total_thr_e;    // Thrown energy [GeV]
-
-  // Loop over event by event
-  for (int ievent = 0; ievent < nentries; ievent++)
-  {
-	t->GetEntry(ievent);
-	
-	Int_t ncluster   = EcalClusters_;
-	Int_t nreconhits = RecoEcalHits_;
-
-	total_cluster_e = 0.0;
-
-	h5->Fill(ncluster, 1.0);
-
-	// Loop over cluster by cluster
-	for (int icluster=0; icluster < ncluster; icluster++)
-	{
-		r = TMath::Sqrt((cluster_x_pos[icluster]*cluster_x_pos[icluster]) + 
-				(cluster_y_pos[icluster]*cluster_y_pos[icluster]) + 
-				(cluster_z_pos[icluster]*cluster_z_pos[icluster]));
-		phi = TMath::ATan(cluster_y_pos[icluster]/cluster_x_pos[icluster]) * TMath::RadToDeg();
-		theta = TMath::ACos(cluster_z_pos[icluster] / r) * TMath::RadToDeg();
-		eta = -1.0 * TMath::Log(TMath::Tan((theta*TMath::DegToRad())/2.0));	
-		cluster_e = cluster_energy[icluster] / 1.e+3;
-		total_cluster_e += cluster_e;
-	}
-	// Writing into ASCII file
-	fprintf(fp,"%d\t%d\t%2.4f\n",ievent,ncluster,total_cluster_e);
-
-	// Select events with one cluster
-        if(ncluster == 1)
-        {       
-		h1->Fill(theta, 1.0);
-                h2->Fill(eta, 1.0);
-                h3->Fill(cluster_e, eta, 1.0);
-                h4->Fill(total_cluster_e, 1.0);
-                h8->Fill(cluster_x_pos[0],cluster_y_pos[0], 1.0);
-
-		if(total_cluster_e > 0.5)
-		{
-			h6->Fill(theta, 1.0);
-			h7->Fill(eta, 1.0);
-		}
-	}
-
-	// Loop over hit by hit
-        for(int ireconhit=0; ireconhit < nreconhits; ireconhit++)
-		h9->Fill(rec_x_pos[ireconhit],rec_y_pos[ireconhit], 1.0);
-  }
-  
-  // Drawing and Saving figures
-  c1->cd();
-  h1->SetLineColor(kBlue);
-  h1->SetLineWidth(2);
-  h1->GetYaxis()->SetRangeUser(0.0,h1->GetMaximum()+50.0);
-  h1->GetXaxis()->SetTitle("#theta [degree]");
-  h1->GetYaxis()->SetTitle("events");
-  h1->GetYaxis()->SetTitleOffset(1.4);
-  h1->DrawClone();
-  c1->SaveAs("results/electron_theta_hist_0GeVto30GeV.png");
-  c1->SaveAs("results/electron_theta_hist_0GeVto30GeV.pdf");
-
-  c2->cd();
-  h2->SetLineColor(kBlue);
-  h2->SetLineWidth(2);
-  h2->GetYaxis()->SetRangeUser(0.0,h2->GetMaximum()+50.0);
-  h2->GetXaxis()->SetTitle("#eta");
-  h2->GetYaxis()->SetTitle("events");
-  h2->GetYaxis()->SetTitleOffset(1.4);
-  h2->DrawClone();
-  c2->SaveAs("results/electron_eta_hist_0GeVto30GeV.png");
-  c2->SaveAs("results/electron_eta_hist_0GeVto30GeV.pdf");
-
-  c3->cd();
-  h3->GetXaxis()->SetTitle("Cluster energy [GeV]");
-  h3->GetYaxis()->SetTitle("#eta");
-  h3->GetYaxis()->SetTitleOffset(1.4);
-  h3->DrawClone("COLZ");
-  c3->SaveAs("results/eletron_E_vs_eta_hist_0GeVto30GeV.png");
-  c3->SaveAs("results/eletron_E_vs_eta_hist_0GeVto30GeV.pdf");
-
-  c4->cd();
-  c4->SetLogy(1);
-  h4->SetLineColor(kBlue);
-  h4->SetLineWidth(2);
-  h4->GetXaxis()->SetTitle("reconstructed energy [GeV]");
-  h4->GetYaxis()->SetTitle("events");
-  h4->GetYaxis()->SetTitleOffset(1.4);
-  h4->DrawClone();
-  c4->SaveAs("results/electron_Erec_hist_0GeVto30GeV.png");
-  c4->SaveAs("results/electron_Erec_hist_0GeVto30GeV.pdf");
-
-  c5->cd();
-  c5->SetLogy(1);
-  h5->SetLineColor(kBlue);
-  h5->SetLineWidth(2);
-  h5->GetXaxis()->SetTitle("Number of Clusters");
-  h5->GetYaxis()->SetTitle("events");
-  h5->GetYaxis()->SetTitleOffset(1.4);
-  h5->DrawClone();
-  c5->SaveAs("results/electron_ncluster_hist_0GeVto30GeV.png");
-  c5->SaveAs("results/electron_ncluster_hist_0GeVto30GeV.pdf");
-
-  c6->cd();
-  h6->SetLineColor(kBlue);
-  h6->SetLineWidth(2);
-  h6->GetYaxis()->SetRangeUser(0.0,h1->GetMaximum()+50.0);
-  h6->GetXaxis()->SetTitle("#theta [degree]");
-  h6->GetYaxis()->SetTitle("events");
-  h6->GetYaxis()->SetTitleOffset(1.4);
-  h6->DrawClone();
-  c6->SaveAs("results/electron_theta_hist_CUT_0GeVto30GeV.png");
-  c6->SaveAs("results/electron_theta_hist_CUT_0GeVto30GeV.pdf");
-
-  c7->cd();
-  h7->SetLineColor(kBlue);
-  h7->SetLineWidth(2);
-  h7->GetYaxis()->SetRangeUser(0.0,h2->GetMaximum()+50.0);
-  h7->GetXaxis()->SetTitle("#eta");
-  h7->GetYaxis()->SetTitle("events");
-  h7->GetYaxis()->SetTitleOffset(1.4);
-  h7->DrawClone();
-  c7->SaveAs("results/electron_eta_hist_CUT_0GeVto30GeV.png");
-  c7->SaveAs("results/electron_eta_hist_CUT_0GeVto30GeV.pdf");
-
-  c8->cd();
-  h8->GetXaxis()->SetTitle("Hit position X [cm]");
-  h8->GetYaxis()->SetTitle("Hit position Y [cm]");
-  h8->GetYaxis()->SetTitleOffset(1.4);
-  h8->DrawClone("COLZ");
-  ell1->Draw("same");
-  ell2->Draw("same");
-  c8->SaveAs("results/electron_hit_pos_cluster_0GeVto30GeV.png");
-  c8->SaveAs("results/electron_hit_pos_cluster_0GeVto30GeV.pdf");
-
-  c9->cd();
-  h9->GetXaxis()->SetTitle("Hit position X [cm]");
-  h9->GetYaxis()->SetTitle("Hit position Y [cm]");
-  h9->GetYaxis()->SetTitleOffset(1.4);
-  h9->DrawClone("COLZ");
-  ell1->Draw("same");
-  ell2->Draw("same");
-  c9->SaveAs("results/electron_hit_pos_all_0GeVto30GeV.png");
-  c9->SaveAs("results/electron_hit_pos_all_0GeVto30GeV.pdf");
-
-  fclose(fp);
-  return 0;
-}
diff --git a/benchmarks/ecal/scripts/deprecated/makeplot_input.C b/benchmarks/ecal/scripts/deprecated/makeplot_input.C
deleted file mode 100644
index 8feab515..00000000
--- a/benchmarks/ecal/scripts/deprecated/makeplot_input.C
+++ /dev/null
@@ -1,102 +0,0 @@
-////////////////////////////////////////
-// Read simulation ROOT input file
-// Write ASCII file
-////////////////////////////////////////
-int makeplot_input(const char* input_fname = "sim_output/sim_emcal_electrons_input.root", const char *fout = "results/sim_1000k_info_1GeV.txt")
-{
-  // Setting figures
-  gROOT->SetStyle("Plain");
-  gStyle->SetLineWidth(3);
-  gStyle->SetOptStat("nem");
-  gStyle->SetPadTickX(1);
-  gStyle->SetPadTickY(1);
-  gStyle->SetPadGridX(1);
-  gStyle->SetPadGridY(1);
-  gStyle->SetPadLeftMargin(0.14);
-  gStyle->SetPadRightMargin(0.14);
-
-  // Output ASCII file
-  FILE *fp = fopen(fout,"a");
-  if(!fp)
-        fprintf(stderr,"error: can't open %s\n", fout);
-
-  // Input ROOT file
-  TFile *f = new TFile(input_fname,"READ");
-  TTree *t = (TTree *)f->Get("events");
-
-  // Set Branch status and addressed
-  t->SetMakeClass(1);
-  t->SetBranchStatus("*", 0);
-
-  Int_t CrystalEcalHits_;
-  t->SetBranchStatus("CrystalEcalHits", 1);
-  t->SetBranchAddress("CrystalEcalHits", &CrystalEcalHits_);
-
-  Int_t mcparticles_;
-  t->SetBranchStatus("mcparticles", 1);
-  t->SetBranchAddress("mcparticles", &mcparticles_);
-
-  const Int_t kMaxCrystalEcalHits = 100000;
-  Double_t truth_deposit[kMaxCrystalEcalHits];
-  Double_t energyDeposit[kMaxCrystalEcalHits];
-  t->SetBranchStatus("CrystalEcalHits.truth.deposit",1);
-  t->SetBranchStatus("CrystalEcalHits.energyDeposit",1);
-  t->SetBranchAddress("CrystalEcalHits.truth.deposit",truth_deposit);
-  t->SetBranchAddress("CrystalEcalHits.energyDeposit",energyDeposit);
-
-  const Int_t kMaxmcparticles = 100000;
-  Double_t px[kMaxmcparticles];
-  Double_t py[kMaxmcparticles];
-  Double_t pz[kMaxmcparticles];
-  Double_t mass[kMaxmcparticles];
-  t->SetBranchStatus("mcparticles.psx",1);
-  t->SetBranchStatus("mcparticles.psy",1);
-  t->SetBranchStatus("mcparticles.psz",1);
-  t->SetBranchStatus("mcparticles.mass",1);
-  t->SetBranchAddress("mcparticles.psx",px);
-  t->SetBranchAddress("mcparticles.psy",py);
-  t->SetBranchAddress("mcparticles.psz",pz);
-  t->SetBranchAddress("mcparticles.mass",mass);
-
-  // Total number of entries
-  Int_t nentries = t->GetEntries();
-
-  // Variables are used in calculation
-  Double_t total_thr_e;       // Thrown energy [GeV]
-  Double_t truth_sim_e;       // truth simulated deposit energy [GeV]
-  Double_t total_truth_sim_e; // Add up truth simulated deposit energy [GeV]
-  Double_t sim_e;             // simulated deposit energy [GeV]
-  Double_t total_sim_e;       // Add up truth simulated deposit energy [GeV]
-  Double_t momentum2;         // momentum squared
-  Double_t mass2;             // mass squared
-
-  // Loop over event by event
-  for (int ievent = 0; ievent < nentries; ievent++)
-  {
-	t->GetEntry(ievent);
-
-        Int_t nCrystalEcalHits = CrystalEcalHits_;
-        Int_t nmcparticle      = 2;
-
-        total_thr_e       = 0.0;
-        total_truth_sim_e = 0.0;
-        total_sim_e       = 0.0;
-
-        momentum2 = px[nmcparticle]*px[nmcparticle]+py[nmcparticle]*py[nmcparticle]+pz[nmcparticle]*pz[nmcparticle];
-        mass2     = mass[nmcparticle]*mass[nmcparticle];
-
-        total_thr_e = sqrt(momentum2 + mass2)/1.e+3;
-
-        // Loop over cluster by cluster
-        for (int isimhit=0; isimhit < nCrystalEcalHits; isimhit++)
-        {
-                truth_sim_e        = truth_deposit[isimhit];
-                sim_e              = energyDeposit[isimhit];
-                total_truth_sim_e += truth_sim_e;
-                total_sim_e       += energyDeposit[isimhit]/1.e+3;
-        }
-        fprintf(fp,"%d\t%2.4f\t%2.4f\t%2.4f\n",ievent,total_thr_e,total_truth_sim_e,total_sim_e);
-  }
-  fclose(fp);
-  return 0;
-}
diff --git a/benchmarks/ecal/scripts/deprecated/makeplot_pi0.C b/benchmarks/ecal/scripts/deprecated/makeplot_pi0.C
deleted file mode 100644
index 94c5bf5b..00000000
--- a/benchmarks/ecal/scripts/deprecated/makeplot_pi0.C
+++ /dev/null
@@ -1,315 +0,0 @@
-////////////////////////////////////////
-// Read reconstruction ROOT output file
-// Plot variables
-////////////////////////////////////////
-int makeplot_pi0(const char* input_fname = "../sim_output/sim_emcal_pi0s_output.root")
-{
-  // Setting figures
-  gROOT->SetStyle("Plain");
-  gStyle->SetLineWidth(3);
-  gStyle->SetOptStat("nem");
-  gStyle->SetPadTickX(1);
-  gStyle->SetPadTickY(1);
-  gStyle->SetPadGridX(1);
-  gStyle->SetPadGridY(1);
-  gStyle->SetPadLeftMargin(0.14);
-  gStyle->SetPadRightMargin(0.14);
-
-  // Input ROOT file
-  TFile *f = new TFile(input_fname,"READ");
-  TTree *t = (TTree *)f->Get("events");
-
-  // Set Branch status and addressed
-  t->SetMakeClass(1);
-  t->SetBranchStatus("*", 0);
-
-  Int_t EcalClusters_;
-  t->SetBranchStatus("EcalClusters", 1);
-  t->SetBranchAddress("EcalClusters", &EcalClusters_);
-
-  Int_t RecoEcalHits_;
-  t->SetBranchStatus("RecoEcalHits", 1);
-  t->SetBranchAddress("RecoEcalHits", &RecoEcalHits_);
-
-  const Int_t kMaxEcalClusters = 100000;
-  Double_t cluster_x_pos[kMaxEcalClusters];
-  Double_t cluster_y_pos[kMaxEcalClusters];
-  Double_t cluster_z_pos[kMaxEcalClusters];
-  Float_t cluster_energy[kMaxEcalClusters];
-  t->SetBranchStatus("EcalClusters.position.x",1);
-  t->SetBranchStatus("EcalClusters.position.y",1);
-  t->SetBranchStatus("EcalClusters.position.z",1);
-  t->SetBranchStatus("EcalClusters.energy",1);
-  t->SetBranchAddress("EcalClusters.position.x",cluster_x_pos);
-  t->SetBranchAddress("EcalClusters.position.y",cluster_y_pos);
-  t->SetBranchAddress("EcalClusters.position.z",cluster_z_pos);
-  t->SetBranchAddress("EcalClusters.energy",cluster_energy);
-
-  const Int_t kMaxRecoEcalHits = 100000;
-  Double_t rec_x_pos[kMaxRecoEcalHits];
-  Double_t rec_y_pos[kMaxRecoEcalHits];
-  Double_t rec_energy[kMaxRecoEcalHits];
-  t->SetBranchStatus("RecoEcalHits.position.x",1);
-  t->SetBranchStatus("RecoEcalHits.position.y",1);
-  t->SetBranchStatus("RecoEcalHits.energy",1);
-  t->SetBranchAddress("RecoEcalHits.position.x",rec_x_pos);
-  t->SetBranchAddress("RecoEcalHits.position.y",rec_y_pos);
-  t->SetBranchAddress("RecoEcalHits.energy",rec_energy);
-
-  // Setting for Canvas
-  TCanvas *c1  = new TCanvas("c1", "c1",  500, 500);
-  TCanvas *c2  = new TCanvas("c2", "c2",  500, 500);
-  TCanvas *c3  = new TCanvas("c3", "c3",  500, 500);
-  TCanvas *c4  = new TCanvas("c4", "c4",  500, 500);
-  TCanvas *c5  = new TCanvas("c5", "c5",  500, 500);
-  TCanvas *c6  = new TCanvas("c6", "c6",  500, 500);
-  TCanvas *c7  = new TCanvas("c7", "c7",  500, 500);
-  TCanvas *c8  = new TCanvas("c8", "c8",  500, 500);
-  TCanvas *c9  = new TCanvas("c9", "c9",  500, 500);
-  TCanvas *c10 = new TCanvas("c10","c10", 500, 500);
-  TCanvas *c11 = new TCanvas("c11","c11", 500, 500);
-  TCanvas *c12 = new TCanvas("c12","c12", 500, 500);
-  TCanvas *c13 = new TCanvas("c13","c13", 500, 500);
-
-  // Declare histograms
-  TH1D *h1  = new TH1D("h1", "Scattering Angle(#theta)",          100,130.0,180.0);
-  TH1D *h2  = new TH1D("h2", "Pseudo-rapidity(#eta)",             100,-5.0,0.0);
-  TH2D *h3  = new TH2D("h3", "Cluster E vs Pseudo-rapidity",      100,-0.5,30.5,100,-5.0,0.0);
-  TH1D *h4  = new TH1D("h4", "Reconstructed energy per event",    100,-0.5,30.5);
-  TH1D *h5  = new TH1D("h5", "Number of Clusters per event",      5,-0.5,4.5);
-  TH1D *h6  = new TH1D("h6", "Scattering Angle(#theta) with CUT", 100,130.0,180.0);
-  TH1D *h7  = new TH1D("h7", "Pseudo-rapidity(#eta) with CUT",    100,-5.0,0.0);
-  TH2D *h8  = new TH2D("h8", "Cluster Hit Position",              62,-62.0,62.0,62,-62.0,62.0);
-  TH2D *h9  = new TH2D("h9", "All Hit Position",                  62,-62.0,62.0,62,-62.0,62.0);
-  TH1D *h10 = new TH1D("h10","Invariant mass",                    60,0.0,300.0);
-  TH1D *h11 = new TH1D("h11","E1",                                100,-0.5,30.5);
-  TH1D *h12 = new TH1D("h12","E2",                                100,-0.5,30.5);
-  TH1D *h13 = new TH1D("h13","angle",                             100,0.0,90.0);
-
-
-  // Declare ellipse for boundary of crystal calorimeter
-  TEllipse *ell1 = new TEllipse(0.0, 0.0, 60.0, 60.0);
-  ell1->SetFillStyle(4000);
-  TEllipse *ell2 = new TEllipse(0.0, 0.0, 12.0, 12.0);
-  ell2->SetFillStyle(4000);
-
-  // Total number of entries
-  Int_t nentries = t->GetEntries();
-
-  // Variables are used in calculation
-  Double_t r;                        // Radius [cm]
-  Double_t phi;                      // Azimuth [degree]
-  Double_t theta;                    // Inclination [degree]
-  Double_t eta;                      // Pseudo-rapidity [unitless]
-  Float_t  cluster_e;                // Cluster energy [GeV]
-  Float_t  total_cluster_e;          // Add up clusters per event [GeV]
-  Double_t dot_product_pos_clusters; // dot product of positions of two photons
-  Double_t mag_pos2_cluster_1;       // squared magnitude of position
-  Double_t mag_pos2_cluster_2;       // squared magnitude of position 
-  Double_t cosine_clusters;          // cos(theta_photons)
-  Double_t theta_photons;            // angle between two photons
-  Double_t invariant_mass;           // M^2 = 2 * p_1 * p_2 * (1 - cos(theta_photons))
-
-  // Loop over event by event
-  for (int ievent = 0; ievent < nentries; ievent++)
-  {
-	t->GetEntry(ievent);
-
-	Int_t ncluster   = EcalClusters_;
-	Int_t nreconhits = RecoEcalHits_;
-
-	total_cluster_e = 0.0;
-
-	h5->Fill(ncluster, 1.0);
-
-	// Loop over cluster by cluster
-	for (int icluster=0; icluster < ncluster; icluster++)
-	{
-		r = TMath::Sqrt((cluster_x_pos[icluster]*cluster_x_pos[icluster]) + 
-				(cluster_y_pos[icluster]*cluster_y_pos[icluster]) + 
-				(cluster_z_pos[icluster]*cluster_z_pos[icluster]));
-		phi = TMath::ATan(cluster_y_pos[icluster]/cluster_x_pos[icluster]) * TMath::RadToDeg();
-		theta = TMath::ACos(cluster_z_pos[icluster] / r) * TMath::RadToDeg();
-		eta = -1.0 * TMath::Log(TMath::Tan((theta*TMath::DegToRad())/2.0));	
-		cluster_e = cluster_energy[icluster] / 1.e+3;
-		total_cluster_e += cluster_e;
-	}
-
-	// Select events with two clusters
-        // To calculate invariant mass
-        // M^2 = 2p1p2(1-cos(theta))
-        // p1 = E1
-        // p2 = E2
-        // theta: angle between two photons
-        if(ncluster == 2)
-        {       
-		h1->Fill(theta, 1.0);
-                h2->Fill(eta, 1.0);
-                h3->Fill(cluster_e, eta, 1.0);
-                h4->Fill(total_cluster_e, 1.0);
-                h8->Fill(cluster_x_pos[0],cluster_y_pos[0], 1.0);
-
-		if(total_cluster_e > 0.5)
-		{
-			h6->Fill(theta, 1.0);
-			h7->Fill(eta, 1.0);
-		}
-
-		// Calculate invariant mass
-                dot_product_pos_clusters = cluster_x_pos[0]*cluster_x_pos[1] + cluster_y_pos[0]*cluster_y_pos[1] + cluster_z_pos[0]*cluster_z_pos[1];
-                mag_pos2_cluster_1 = (cluster_x_pos[0]*cluster_x_pos[0]) + (cluster_y_pos[0]*cluster_y_pos[0]) + (cluster_z_pos[0]*cluster_z_pos[0]);
-                mag_pos2_cluster_2 = (cluster_x_pos[1]*cluster_x_pos[1]) + (cluster_y_pos[1]*cluster_y_pos[1]) + (cluster_z_pos[1]*cluster_z_pos[1]);
-                cosine_clusters = (dot_product_pos_clusters/TMath::Sqrt(mag_pos2_cluster_1*mag_pos2_cluster_2));
-                theta_photons = TMath::ACos(cosine_clusters)*TMath::RadToDeg();
-                
-                invariant_mass = TMath::Sqrt(2.0*cluster_energy[0]*cluster_energy[1]*(1.0 - cosine_clusters));
-                
-                // Fill histograms
-                h10->Fill(invariant_mass, 1.0);
-                h11->Fill(cluster_energy[0], 1.0);
-                h12->Fill(cluster_energy[1], 1.0);
-                h13->Fill(theta_photons, 1.0);
-	}
-
-	// Loop over hit by hit
-        for(int ireconhit=0; ireconhit < nreconhits; ireconhit++)
-		h9->Fill(rec_x_pos[ireconhit],rec_y_pos[ireconhit], 1.0);
-  }
-
-  // Drawing and Saving figures
-  c1->cd();
-  h1->SetLineColor(kBlue);
-  h1->SetLineWidth(2);
-  h1->GetYaxis()->SetRangeUser(0.0,h1->GetMaximum()+10.0);
-  h1->GetXaxis()->SetTitle("#theta [degree]");
-  h1->GetYaxis()->SetTitle("events");
-  h1->GetYaxis()->SetTitleOffset(1.4);
-  h1->DrawClone();
-  c1->SaveAs("results/pi0_theta_hist_0GeVto30GeV.png");
-  c1->SaveAs("results/pi0_theta_hist_0GeVto30GeV.pdf");
-
-  c2->cd();
-  h2->SetLineColor(kBlue);
-  h2->SetLineWidth(2);
-  h2->GetYaxis()->SetRangeUser(0.0,h2->GetMaximum()+10.0);
-  h2->GetXaxis()->SetTitle("#eta");
-  h2->GetYaxis()->SetTitle("events");
-  h2->GetYaxis()->SetTitleOffset(1.4);
-  h2->DrawClone();
-  c2->SaveAs("results/pi0_eta_hist_0GeVto30GeV.png");
-  c2->SaveAs("results/pi0_eta_hist_0GeVto30GeV.pdf");
-
-  c3->cd();
-  h3->GetXaxis()->SetTitle("Cluster energy [GeV]");
-  h3->GetYaxis()->SetTitle("#eta");
-  h3->GetYaxis()->SetTitleOffset(1.4);
-  h3->DrawClone("COLZ");
-  c3->SaveAs("results/pi0_E_vs_eta_hist_0GeVto30GeV.png");
-  c3->SaveAs("results/pi0_E_vs_eta_hist_0GeVto30GeV.pdf");
-
-  c4->cd();
-  c4->SetLogy(1);
-  h4->SetLineColor(kBlue);
-  h4->SetLineWidth(2);
-  h4->GetXaxis()->SetTitle("reconstructed energy [GeV]");
-  h4->GetYaxis()->SetTitle("events");
-  h4->GetYaxis()->SetTitleOffset(1.4);
-  h4->DrawClone();
-  c4->SaveAs("results/pi0_Erec_hist_0GeVto30GeV.png");
-  c4->SaveAs("results/pi0_Erec_hist_0GeVto30GeV.pdf");
-
-  c5->cd();
-  c5->SetLogy(1);
-  h5->SetLineColor(kBlue);
-  h5->SetLineWidth(2);
-  h5->GetXaxis()->SetTitle("Number of Clusters");
-  h5->GetYaxis()->SetTitle("events");
-  h5->GetYaxis()->SetTitleOffset(1.4);
-  h5->DrawClone();
-  c5->SaveAs("results/pi0_ncluster_hist_0GeVto30GeV.png");
-  c5->SaveAs("results/pi0_ncluster_hist_0GeVto30GeV.pdf");
-
-  c6->cd();
-  h6->SetLineColor(kBlue);
-  h6->SetLineWidth(2);
-  h6->GetYaxis()->SetRangeUser(0.0,h1->GetMaximum()+10.0);
-  h6->GetXaxis()->SetTitle("#theta [degree]");
-  h6->GetYaxis()->SetTitle("events");
-  h6->GetYaxis()->SetTitleOffset(1.4);
-  h6->DrawClone();
-  c6->SaveAs("results/pi0_theta_hist_CUT_0GeVto30GeV.png");
-  c6->SaveAs("results/pi0_theta_hist_CUT_0GeVto30GeV.pdf");
-
-  c7->cd();
-  h7->SetLineColor(kBlue);
-  h7->SetLineWidth(2);
-  h7->GetYaxis()->SetRangeUser(0.0,h2->GetMaximum()+10.0);
-  h7->GetXaxis()->SetTitle("#eta");
-  h7->GetYaxis()->SetTitle("events");
-  h7->GetYaxis()->SetTitleOffset(1.4);
-  h7->DrawClone();
-  c7->SaveAs("results/pi0_eta_hist_CUT_0GeVto30GeV.png");
-  c7->SaveAs("results/pi0_eta_hist_CUT_0GeVto30GeV.pdf");
-
-  c8->cd();
-  h8->GetXaxis()->SetTitle("Hit position X [cm]");
-  h8->GetYaxis()->SetTitle("Hit position Y [cm]");
-  h8->GetYaxis()->SetTitleOffset(1.4);
-  h8->DrawClone("COLZ");
-  ell1->Draw("same");
-  ell2->Draw("same");
-  c8->SaveAs("results/pi0_hit_pos_cluster_0GeVto30GeV.png");
-  c8->SaveAs("results/pi0_hit_pos_cluster_0GeVto30GeV.pdf");
-
-  c9->cd();
-  h9->GetXaxis()->SetTitle("Hit position X [cm]");
-  h9->GetYaxis()->SetTitle("Hit position Y [cm]");
-  h9->GetYaxis()->SetTitleOffset(1.4);
-  h9->DrawClone("COLZ");
-  ell1->Draw("same");
-  ell2->Draw("same");
-  c9->SaveAs("results/pi0_hit_pos_all_0GeVto30GeV.png");
-  c9->SaveAs("results/pi0_hit_pos_all_0GeVto30GeV.pdf");
-
-  c10->cd();
-  h10->SetLineColor(kBlue);
-  h10->SetLineWidth(2);
-  h10->GetXaxis()->SetTitle("Invariant mass [MeV]");
-  h10->GetYaxis()->SetTitle("events");
-  h10->GetYaxis()->SetTitleOffset(1.4);
-  h10->DrawClone();
-  c10->SaveAs("results/pi0_invariant_mass_hist_0GeVto30GeV.png"); 
-  c10->SaveAs("results/pi0_invariant_mass_hist_0GeVto30GeV.pdf");
-  
-  c11->cd();
-  h11->SetLineColor(kBlue);
-  h11->SetLineWidth(2);
-  h11->GetXaxis()->SetTitle("Cluster energy 1 [MeV]");
-  h11->GetYaxis()->SetTitle("events");
-  h11->GetYaxis()->SetTitleOffset(1.4);
-  h11->DrawClone();
-  c11->SaveAs("results/pi0_E1_hist_0GeVto30GeV.png");
-  c11->SaveAs("results/pi0_E1_hist_0GeVto30GeV.pdf");
-  
-  c12->cd();
-  h12->SetLineColor(kBlue);
-  h12->SetLineWidth(2);
-  h12->GetXaxis()->SetTitle("Cluster energy 2 [MeV]");
-  h12->GetYaxis()->SetTitle("events");
-  h12->GetYaxis()->SetTitleOffset(1.4);
-  h12->DrawClone();
-  c12->SaveAs("results/pi0_E2_hist_0GeVto30GeV.png");
-  c12->SaveAs("results/pi0_E2_hist_0GeVto30GeV.pdf");
-
-  c13->cd();
-  h13->SetLineColor(kBlue);
-  h13->SetLineWidth(2);
-  h13->GetXaxis()->SetTitle("Angle between two photons [degree]");
-  h13->GetYaxis()->SetTitle("events");
-  h13->GetYaxis()->SetTitleOffset(1.4);
-  h13->DrawClone();
-  c13->SaveAs("results/pi0_angle_2photons_0GeVto30GeV.png");
-  c13->SaveAs("results/pi0_angle_2photons_0GeVto30GeV.pdf");
-
-  return 0;
-}
diff --git a/benchmarks/ecal/scripts/deprecated/rdf_test.cxx b/benchmarks/ecal/scripts/deprecated/rdf_test.cxx
deleted file mode 100644
index 7dc3b14e..00000000
--- a/benchmarks/ecal/scripts/deprecated/rdf_test.cxx
+++ /dev/null
@@ -1,66 +0,0 @@
-#include "ROOT/RDataFrame.hxx"
-#include <iostream>
-
-#include "dd4pod/Geant4ParticleCollection.h"
-
-//namespace edm4hep {
-//
-//std::vector<float> pt (std::vector<MCParticleData> const& in){
-//  std::vector<float> result;
-//  for (size_t i = 0; i < in.size(); ++i) {
-//    result.push_back(std::sqrt(in[i].momentum.x * in[i].momentum.x + in[i].momentum.y * in[i].momentum.y));
-//  }
-//  return result;
-//}
-//
-//std::vector<float> eta(std::vector<MCParticleData> const& in){
-//  std::vector<float> result;
-//  ROOT::Math::PxPyPzMVector lv;
-//  for (size_t i = 0; i < in.size(); ++i) {
-//    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
-//    result.push_back(lv.Eta());
-//  }
-//  return result;
-//}
-//
-//std::vector<float> cos_theta(std::vector<MCParticleData> const& in){
-//  std::vector<float> result;
-//  ROOT::Math::PxPyPzMVector lv;
-//  for (size_t i = 0; i < in.size(); ++i) {
-//    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
-//    result.push_back(cos(lv.Theta()));
-//  }
-//  return result;
-//}
-//
-//}
-
-
-std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
-  std::vector<float> result;
-  for (size_t i = 0; i < in.size(); ++i) {
-    result.push_back(std::sqrt(in[i].psx * in[i].psx + in[i].psy * in[i].psy));
-  }
-  return result;
-}
-
-int rdf_test() {
-
-  ROOT::EnableImplicitMT();
-
-  ROOT::RDataFrame df("events", "sim_barrel_clusters.root");
-
-  auto df2 =
-      df.Define("MCParticles_pt", pt, {"mcparticles"});
-          //.Define("MCParticles_eta", edm4hep::eta, {"mcparticles"})
-          //.Define("MCParticles_cosTheta", edm4hep::cos_theta, {"mcparticles"});
-
-  std::string outfilename = "rdf_test.root";
-
-  df2.Snapshot("events", outfilename,
-               {
-                   "MCParticles_pt",
-                   "mcparticles"
-               });
-  return 0;
-}
diff --git a/benchmarks/ecal/scripts/deprecated/read_eng.C b/benchmarks/ecal/scripts/deprecated/read_eng.C
deleted file mode 100644
index cedc16da..00000000
--- a/benchmarks/ecal/scripts/deprecated/read_eng.C
+++ /dev/null
@@ -1,70 +0,0 @@
-int read_eng(const char* output_fname = "sim_output/read_eng_output.root", const char *fname = "results/eng.txt")
-{
-  // Input ASCII File
-  FILE *fin = fopen(fname,"r");
-
-  if(!fin)
-	fprintf(stderr,"error: can't open %s\n", fname);
-
-  // Output ROOT File
-  TFile *rootFile = new TFile(output_fname,"recreate");
-  // ROOT tree that goes into this ROOT file
-  TTree *tEngRes = new TTree("tEngRes","Energy Resolution Information");
-  // Variables that go into ROOT file
-  Int_t    ievent;
-  Double_t tot_clust_eng; // [GeV]
-  Int_t    ncluster;
-  Double_t mc_eng;        // [GeV]
-  Double_t sim_tru_eng;   // [GeV]
-  Double_t sim_eng;       // [GeV]
-  // Declare ROOT tree branches
-  tEngRes->Branch("ievent",        &ievent,        "ievent/I");
-  tEngRes->Branch("tot_clust_eng", &tot_clust_eng, "tot_clust_eng/D");
-  tEngRes->Branch("ncluster",      &ncluster,      "ncluster/I");
-  tEngRes->Branch("mc_eng",        &mc_eng,        "mc_eng/D");
-  tEngRes->Branch("sim_tru_eng",   &sim_tru_eng,   "sim_tru_eng/D");
-  tEngRes->Branch("sim_eng",       &sim_eng,       "sim_eng/D");
-
-  char inbuf[1024];
-  int  nevents_readout = 0;
-  int  linenum = 0;
-
-  TString s;
-  TObjArray *o = 0;
-
-  // Read each line
-  while(fgets(inbuf,1024,fin))
-  {
-	nevents_readout++;
-	s = inbuf;
-	o = s.Tokenize(" \t\n");
-	for(int icol = 0; icol < o->GetEntries(); icol++)
-	{
-		s = ((TObjString *)(*o)[icol])->GetString();
-		if(icol==0)
-			sscanf(s.Data(),"%d", &ievent);
-		if(icol==1)
-			sscanf(s.Data(),"%lf", &mc_eng);
-		if(icol==2)
-			sscanf(s.Data(),"%lf", &sim_tru_eng);
-		if(icol==3)
-			sscanf(s.Data(),"%lf", &sim_eng);
-		if(icol==5)
-			sscanf(s.Data(),"%d", &ncluster);
-		if(icol==6)
-			sscanf(s.Data(),"%lf", &tot_clust_eng);
-	}
-	tEngRes->Fill();
-	// Clean up
-	delete o;
-  } 
-  // Done reading ASCII file. So finialize the ROOT tree and close the ROOT file
-  fprintf(stdout,"%d events read in!\n",nevents_readout);
-  // Finalize ROOT file
-  tEngRes->Write();
-  // Close ROOT file
-  rootFile->Close();
-  fclose(fin);
-
-  return 0;
-}
diff --git a/benchmarks/ecal/scripts/deprecated/rec_emcal_electrons_reader.C b/benchmarks/ecal/scripts/deprecated/rec_emcal_electrons_reader.C
deleted file mode 100644
index 38d59ba3..00000000
--- a/benchmarks/ecal/scripts/deprecated/rec_emcal_electrons_reader.C
+++ /dev/null
@@ -1,370 +0,0 @@
-////////////////////////////////////////
-// Read reconstruction ROOT output file
-// Plot variables
-////////////////////////////////////////
-int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const char* input_fname = "sim_output/sim_emcal_electrons_output.root")
-{
-  // Setting figures
-  gROOT->SetStyle("Plain");
-  gStyle->SetLineWidth(3);
-  gStyle->SetOptStat("nem");
-  gStyle->SetPadTickX(1);
-  gStyle->SetPadTickY(1);
-  gStyle->SetPadGridX(1);
-  gStyle->SetPadGridY(1);
-  gStyle->SetPadLeftMargin(0.14);
-  gStyle->SetPadRightMargin(0.14);
-
-  // Input ROOT file
-  TFile *f = new TFile(input_fname,"READ");
-  TTree *t = (TTree *)f->Get("events");
-
-  // Set Branch status and addressed
-  t->SetMakeClass(1);
-  t->SetBranchStatus("*", 0);
-/*
-  Int_t mcparticles2_;
-  t->SetBranchStatus("mcparticles2", 1);
-  t->SetBranchAddress("mcparticles2", &mcparticles2_);
-  Int_t CrystalEcalHits2_;
-  t->SetBranchStatus("CrystalEcalHits2", 1);
-  t->SetBranchAddress("CrystalEcalHits2", &CrystalEcalHits2_);
-*/
-  Int_t RecoEcalHits_;
-  t->SetBranchStatus("RecoEcalHits", 1);
-  t->SetBranchAddress("RecoEcalHits", &RecoEcalHits_);
-  Int_t EcalClusters_;
-  t->SetBranchStatus("EcalClusters", 1);
-  t->SetBranchAddress("EcalClusters", &EcalClusters_);
-/*
-  const Int_t kMaxmcparticles2 = 100000;
-  Double_t px[kMaxmcparticles2];
-  Double_t py[kMaxmcparticles2];
-  Double_t pz[kMaxmcparticles2];
-  Double_t mass[kMaxmcparticles2];
-  t->SetBranchStatus("mcparticles2.psx",1);
-  t->SetBranchStatus("mcparticles2.psy",1);
-  t->SetBranchStatus("mcparticles2.psz",1);
-  t->SetBranchStatus("mcparticles2.mass",1);
-  t->SetBranchAddress("mcparticles2.psx",px);
-  t->SetBranchAddress("mcparticles2.psy",py);
-  t->SetBranchAddress("mcparticles2.psz",pz);
-  t->SetBranchAddress("mcparticles2.mass",mass);
-  const Int_t kMaxCrystalEcalHits2 = 100000;
-  Double_t truth_deposit[kMaxCrystalEcalHits2];
-  Double_t energyDeposit[kMaxCrystalEcalHits2];
-  t->SetBranchStatus("CrystalEcalHits2.truth.deposit",1);
-  t->SetBranchStatus("CrystalEcalHits2.energyDeposit",1);
-  t->SetBranchAddress("CrystalEcalHits2.truth.deposit",truth_deposit);
-  t->SetBranchAddress("CrystalEcalHits2.energyDeposit",energyDeposit);
-*/
-  const Int_t kMaxRecoEcalHits = 100000;
-  Double_t rec_x_pos[kMaxRecoEcalHits];
-  Double_t rec_y_pos[kMaxRecoEcalHits];
-  Double_t rec_energy[kMaxRecoEcalHits];
-  t->SetBranchStatus("RecoEcalHits.position.x",1);
-  t->SetBranchStatus("RecoEcalHits.position.y",1);
-  t->SetBranchStatus("RecoEcalHits.energy",1);
-  t->SetBranchAddress("RecoEcalHits.position.x",rec_x_pos);
-  t->SetBranchAddress("RecoEcalHits.position.y",rec_y_pos);
-  t->SetBranchAddress("RecoEcalHits.energy",rec_energy);
-  const Int_t kMaxEcalClusters = 100000;
-  Double_t cluster_x_pos[kMaxEcalClusters];
-  Double_t cluster_y_pos[kMaxEcalClusters];
-  Double_t cluster_z_pos[kMaxEcalClusters];
-  Float_t cluster_energy[kMaxEcalClusters];
-  t->SetBranchStatus("EcalClusters.position.x",1);
-  t->SetBranchStatus("EcalClusters.position.y",1);
-  t->SetBranchStatus("EcalClusters.position.z",1);
-  t->SetBranchStatus("EcalClusters.energy",1);
-  t->SetBranchAddress("EcalClusters.position.x",cluster_x_pos);
-  t->SetBranchAddress("EcalClusters.position.y",cluster_y_pos);
-  t->SetBranchAddress("EcalClusters.position.z",cluster_z_pos);
-  t->SetBranchAddress("EcalClusters.energy",cluster_energy);
-
-  // Setting for Canvas
-  TCanvas *c1  = new TCanvas("c1",  "c1",  500, 500);
-  TCanvas *c2  = new TCanvas("c2",  "c2",  500, 500);
-  TCanvas *c3  = new TCanvas("c3",  "c3",  500, 500);
-  TCanvas *c4  = new TCanvas("c4",  "c4",  500, 500);
-  TCanvas *c5  = new TCanvas("c5",  "c5",  500, 500);
-  TCanvas *c6  = new TCanvas("c6",  "c6",  500, 500);
-  TCanvas *c7  = new TCanvas("c7",  "c7",  500, 500);
-  TCanvas *c8  = new TCanvas("c8",  "c8",  500, 500);
-  TCanvas *c9  = new TCanvas("c9",  "c9",  500, 500);
-  TCanvas *c10 = new TCanvas("c10", "c10", 700, 500);
-  TCanvas *c11 = new TCanvas("c11", "c11", 700, 500);
-  TCanvas *c12 = new TCanvas("c12", "c12", 700, 500);
-  TCanvas *c13 = new TCanvas("c13", "c13", 700, 500);
-  TCanvas *c14 = new TCanvas("c14", "c14", 700, 500);
-
-  // Declare histograms
-  TH1D *h1  = new TH1D("h1","Scattering Angle(#theta)",                    100,130.0,180.0);
-  TH1D *h2  = new TH1D("h2","Pseudo-rapidity(#eta)",                       100,-5.0,0.0);
-  TH2D *h3  = new TH2D("h3","Cluster E vs Pseudo-rapidity",                100,e_start-0.5,e_end+0.5,100,-5.0,0.0);
-  TH1D *h4  = new TH1D("h4","Reconstructed energy per event",              100,e_start-0.5,e_end+0.5);
-  TH1D *h5  = new TH1D("h5","Number of Clusters per event",                5,-0.5,4.5);
-  TH1D *h6  = new TH1D("h6","Scattering Angle(#theta) with CUT",           100,130.0,180.0);
-  TH1D *h7  = new TH1D("h7","Pseudo-rapidity(#eta) with CUT",              100,-5.0,0.0);
-  TH2D *h8  = new TH2D("h8","Cluster Hit Position",                        62,-62.0,62.0,62,-62.0,62.0);
-  TH2D *h9  = new TH2D("h9","All Hit Position",                            62,-62.0,62.0,62,-62.0,62.0);
-  TH1D *h10 = new TH1D("hEnergyRes","Energy Resolution",                   100,-0.3,0.3);
-  TH1D *h11 = new TH1D("hEnergyResCUT","Energy Resolution with CUT",       100,-0.3,0.3);
-  TH1D *h12 = new TH1D("h12","Thrown momentum",                            61,e_start-0.5,e_end+0.5);
-  TH1D *h13 = new TH1D("h13","Thrown momentum for reconstructed particle", 61,e_start-0.5,e_end+0.5);
-  TH1D *h14 = new TH1D("h14","Ratio p_{rec}/p_{thr}",                      61,e_start-0.5,e_end+0.5);
-
-  // Declare ellipse for boundary of crystal calorimeter
-  TEllipse *ell1 = new TEllipse(0.0, 0.0, 60.0, 60.0);
-  ell1->SetFillStyle(4000);
-  TEllipse *ell2 = new TEllipse(0.0, 0.0, 12.0, 12.0);
-  ell2->SetFillStyle(4000);
-
-  // Total number of entries
-  Int_t nentries = t->GetEntries();
-
-  // Variables are used in calculation
-  Double_t r;                 // Radius [cm]
-  Double_t phi;               // Azimuth [degree]
-  Double_t theta;             // Inclination [degree]
-  Double_t eta;               // Pseudo-rapidity [unitless]
-  Float_t  cluster_e;         // Cluster energy [GeV]
-  Float_t  total_cluster_e;   // Add up clusters per event [GeV]
-  Double_t total_thr_e;       // Thrown energy [GeV]
-  Double_t total_truth_sim_e; // Add up truth simulated deposit energy [GeV]
-  Double_t total_sim_e;       // Add up truth simulated deposit energy [GeV]
-  Double_t momentum2;         // momentum squared
-  Double_t momentum;          // momentum
-  Double_t mass2;             // mass squared
-  Double_t diff_energy;       // delta(E) = reconstructed energy - thrown energy 
-  Double_t eng_res;           // delta(E)/thrown energy
-
-  // Loop over event by event
-  for (int ievent = 0; ievent < nentries; ievent++)
-  {
-	// Read event by event
-	t->GetEntry(ievent);
-	// Read number of hits/clusters
-/*
-	Int_t nmcparticle      = 2;
-	Int_t nCrystalEcalHits = CrystalEcalHits2_;
-*/
-	Int_t nreconhits       = RecoEcalHits_;
-	Int_t ncluster         = EcalClusters_;
-	// Initialize total energy variables
-	total_thr_e       = 0.0;
-        total_truth_sim_e = 0.0;
-        total_sim_e       = 0.0;
-	total_cluster_e   = 0.0;
-	// Thrown energy, momentum, and mass
-/*
-	momentum2   = px[nmcparticle]*px[nmcparticle]+py[nmcparticle]*py[nmcparticle]+pz[nmcparticle]*pz[nmcparticle];
-        momentum    = TMath::Sqrt(momentum2);
-	mass2       = mass[nmcparticle]*mass[nmcparticle];
-        total_thr_e = sqrt(momentum2 + mass2)/1.e+3;
-	h12->Fill(momentum, 1.0);
-
-        // Loop over simulated hit by simulated hit
-        for (int isimhit=0; isimhit < nCrystalEcalHits; isimhit++)
-        {
-                total_truth_sim_e += truth_deposit[isimhit];
-                total_sim_e       += energyDeposit[isimhit]/1.e+3;
-        }
-*/        
-	// Loop over reconstructed hit by reconstructed hit
-        for(int ireconhit=0; ireconhit < nreconhits; ireconhit++)
-                h9->Fill(rec_x_pos[ireconhit],rec_y_pos[ireconhit], 1.0);
-
-	h5->Fill(ncluster, 1.0);
-
-	// Loop over cluster by cluster
-	for (int icluster=0; icluster < ncluster; icluster++)
-	{
-		r = TMath::Sqrt((cluster_x_pos[icluster]*cluster_x_pos[icluster]) + 
-				(cluster_y_pos[icluster]*cluster_y_pos[icluster]) + 
-				(cluster_z_pos[icluster]*cluster_z_pos[icluster]));
-		phi = TMath::ATan(cluster_y_pos[icluster]/cluster_x_pos[icluster]) * TMath::RadToDeg();
-		theta = TMath::ACos(cluster_z_pos[icluster] / r) * TMath::RadToDeg();
-		eta = -1.0 * TMath::Log(TMath::Tan((theta*TMath::DegToRad())/2.0));	
-		cluster_e = cluster_energy[icluster] / 1.e+3;
-		total_cluster_e += cluster_e;
-	}
-
-	// Select events with one cluster
-        if(ncluster == 1)
-        {       
-                diff_energy = total_cluster_e - total_thr_e;
-                eng_res     = (diff_energy / total_thr_e);
-                h10->Fill(eng_res, 1.0);
-
-		h1->Fill(theta, 1.0);
-                h2->Fill(eta, 1.0);
-                h3->Fill(cluster_e, eta, 1.0);
-                h4->Fill(total_cluster_e, 1.0);
-                h8->Fill(cluster_x_pos[0],cluster_y_pos[0], 1.0);
-		h10->Fill(eng_res, 1.0);
-
-		if(total_cluster_e > 0.2*total_thr_e)
-		{
-			h6->Fill(theta, 1.0);
-			h7->Fill(eta, 1.0);
-			h11->Fill(eng_res, 1.0);
-		}
-/*
-		if(total_cluster_e > 0.9*total_thr_e)
-			h13->Fill(momentum, 1.0);
-*/
-	}
-  } 
-  // Drawing and Saving figures
-  c1->cd();
-  h1->SetLineColor(kBlue);
-  h1->SetLineWidth(2);
-  h1->GetYaxis()->SetRangeUser(0.0,h1->GetMaximum()+50.0);
-  h1->GetXaxis()->SetTitle("#theta [degree]");
-  h1->GetYaxis()->SetTitle("events");
-  h1->GetYaxis()->SetTitleOffset(1.4);
-  h1->DrawClone();
-  c1->SaveAs("results/electron_theta_hist_0GeVto30GeV.png");
-  c1->SaveAs("results/electron_theta_hist_0GeVto30GeV.pdf");
-
-  c2->cd();
-  h2->SetLineColor(kBlue);
-  h2->SetLineWidth(2);
-  h2->GetYaxis()->SetRangeUser(0.0,h2->GetMaximum()+50.0);
-  h2->GetXaxis()->SetTitle("#eta");
-  h2->GetYaxis()->SetTitle("events");
-  h2->GetYaxis()->SetTitleOffset(1.4);
-  h2->DrawClone();
-  c2->SaveAs("results/electron_eta_hist_0GeVto30GeV.png");
-  c2->SaveAs("results/electron_eta_hist_0GeVto30GeV.pdf");
-
-  c3->cd();
-  h3->GetXaxis()->SetTitle("Cluster energy [GeV]");
-  h3->GetYaxis()->SetTitle("#eta");
-  h3->GetYaxis()->SetTitleOffset(1.4);
-  h3->DrawClone("COLZ");
-  c3->SaveAs("results/eletron_E_vs_eta_hist_0GeVto30GeV.png");
-  c3->SaveAs("results/eletron_E_vs_eta_hist_0GeVto30GeV.pdf");
-
-  c4->cd();
-  c4->SetLogy(1);
-  h4->SetLineColor(kBlue);
-  h4->SetLineWidth(2);
-  h4->GetXaxis()->SetTitle("reconstructed energy [GeV]");
-  h4->GetYaxis()->SetTitle("events");
-  h4->GetYaxis()->SetTitleOffset(1.4);
-  h4->DrawClone();
-  c4->SaveAs("results/electron_Erec_hist_0GeVto30GeV.png");
-  c4->SaveAs("results/electron_Erec_hist_0GeVto30GeV.pdf");
-
-  c5->cd();
-  c5->SetLogy(1);
-  h5->SetLineColor(kBlue);
-  h5->SetLineWidth(2);
-  h5->GetXaxis()->SetTitle("Number of Clusters");
-  h5->GetYaxis()->SetTitle("events");
-  h5->GetYaxis()->SetTitleOffset(1.4);
-  h5->DrawClone();
-  c5->SaveAs("results/electron_ncluster_hist_0GeVto30GeV.png");
-  c5->SaveAs("results/electron_ncluster_hist_0GeVto30GeV.pdf");
-
-  c6->cd();
-  h6->SetLineColor(kBlue);
-  h6->SetLineWidth(2);
-  h6->GetYaxis()->SetRangeUser(0.0,h1->GetMaximum()+50.0);
-  h6->GetXaxis()->SetTitle("#theta [degree]");
-  h6->GetYaxis()->SetTitle("events");
-  h6->GetYaxis()->SetTitleOffset(1.4);
-  h6->DrawClone();
-  c6->SaveAs("results/electron_theta_hist_CUT_0GeVto30GeV.png");
-  c6->SaveAs("results/electron_theta_hist_CUT_0GeVto30GeV.pdf");
-
-  c7->cd();
-  h7->SetLineColor(kBlue);
-  h7->SetLineWidth(2);
-  h7->GetYaxis()->SetRangeUser(0.0,h2->GetMaximum()+50.0);
-  h7->GetXaxis()->SetTitle("#eta");
-  h7->GetYaxis()->SetTitle("events");
-  h7->GetYaxis()->SetTitleOffset(1.4);
-  h7->DrawClone();
-  c7->SaveAs("results/electron_eta_hist_CUT_0GeVto30GeV.png");
-  c7->SaveAs("results/electron_eta_hist_CUT_0GeVto30GeV.pdf");
-
-  c8->cd();
-  h8->GetXaxis()->SetTitle("Hit position X [cm]");
-  h8->GetYaxis()->SetTitle("Hit position Y [cm]");
-  h8->GetYaxis()->SetTitleOffset(1.4);
-  h8->DrawClone("COLZ");
-  ell1->Draw("same");
-  ell2->Draw("same");
-  c8->SaveAs("results/electron_hit_pos_cluster_0GeVto30GeV.png");
-  c8->SaveAs("results/electron_hit_pos_cluster_0GeVto30GeV.pdf");
-
-  c9->cd();
-  h9->GetXaxis()->SetTitle("Hit position X [cm]");
-  h9->GetYaxis()->SetTitle("Hit position Y [cm]");
-  h9->GetYaxis()->SetTitleOffset(1.4);
-  h9->DrawClone("COLZ");
-  ell1->Draw("same");
-  ell2->Draw("same");
-  c9->SaveAs("results/electron_hit_pos_all_0GeVto30GeV.png");
-  c9->SaveAs("results/electron_hit_pos_all_0GeVto30GeV.pdf");
-
-  c10->cd();
-  h10->SetLineColor(kBlue);
-  h10->SetLineWidth(2);
-  h10->GetXaxis()->SetTitle("#DeltaE/E");
-  h10->GetYaxis()->SetTitle("events");
-  h10->GetYaxis()->SetTitleOffset(1.4);
-  h10->Fit("gaus");
-  h10->DrawClone();
-  c10->SaveAs("results/hEngRes.png");
-  c10->SaveAs("results/hEngRes.pdf");
-
-  c11->cd();
-  h11->SetLineColor(kBlue);
-  h11->SetLineWidth(2);
-  h11->GetXaxis()->SetTitle("#DeltaE/E");
-  h11->GetYaxis()->SetTitle("events");
-  h11->GetYaxis()->SetTitleOffset(1.4);
-  h11->Fit("gaus");
-  h11->DrawClone();
-  c11->SaveAs("results/hEngRes_CUT.png");
-  c11->SaveAs("results/hEngRes_CUT.pdf");
-
-  c12->cd();
-  h12->SetLineColor(kBlue);
-  h12->SetLineWidth(2);
-  h12->GetXaxis()->SetTitle("Reconstructed Momentum [GeV]");
-  h12->GetYaxis()->SetTitle("events");
-  h12->GetYaxis()->SetTitleOffset(1.4);
-  h12->DrawClone();
-  c12->SaveAs("results/hThrMom.png");
-  c12->SaveAs("results/hThrMom.pdf");
-
-  c13->cd();
-  h13->SetLineColor(kBlue);
-  h13->SetLineWidth(2);
-  h13->GetXaxis()->SetTitle("Thrown Momentum [GeV]");
-  h13->GetYaxis()->SetTitle("events");
-  h13->GetYaxis()->SetTitleOffset(1.4);
-  h13->DrawClone();
-  c13->SaveAs("results/hThrMom_Rec.png");
-  c13->SaveAs("results/hThrMom_Rec.pdf");
-
-  c14->cd();
-  h14 = (TH1D*)h12->Clone();
-  h14->Divide(h13);
-  h14->SetTitle("Ratio p_{rec}/p_{thr}");
-  h14->SetLineColor(kBlue);
-  h14->SetLineWidth(2);
-  h14->GetXaxis()->SetTitle("p [GeV]");
-  h14->GetYaxis()->SetTitle("p_{rec}/p_{thr}");
-  h14->GetYaxis()->SetTitleOffset(1.4);
-  h14->DrawClone();
-  c14->SaveAs("results/hAcceptance.png");
-  c14->SaveAs("results/hAcceptance.pdf");
-  
-  return 0;
-}
diff --git a/benchmarks/ecal/scripts/deprecated/rec_emcal_resolution_analysis.cxx b/benchmarks/ecal/scripts/deprecated/rec_emcal_resolution_analysis.cxx
deleted file mode 100644
index 7fe481d8..00000000
--- a/benchmarks/ecal/scripts/deprecated/rec_emcal_resolution_analysis.cxx
+++ /dev/null
@@ -1,248 +0,0 @@
-////////////////////////////////////////
-// Read reconstruction ROOT output file
-// Plot energy resolution
-// 2021/01/20
-////////////////////////////////////////
-
-#include "ROOT/RDataFrame.hxx"
-#include <iostream>
-
-#include "dd4pod/CalorimeterHitCollection.h"
-#include "dd4pod/Geant4ParticleCollection.h"
-#include "dd4pod/TrackerHitCollection.h"
-#include "eicd/CalorimeterHitCollection.h"
-#include "eicd/CalorimeterHitData.h"
-#include "eicd/ClusterCollection.h"
-#include "eicd/ClusterData.h"
-
-#include "TCanvas.h"
-#include "TH1.h"
-#include "TH1D.h"
-#include "TMath.h"
-#include "TStyle.h"
-
-// If root cannot find a dd4pod header make sure to add the include dir to ROOT_INCLUDE_PATH:
-//  export ROOT_INCLUDE_PATH=$HOME/include:$ROOT_INCLUDE_PATH
-//  export ROOT_INCLUDE_PATH=/home/jihee/stow/development/include:$ROOT_INCLUDE_PATH
-
-using ROOT::RDataFrame;
-using namespace ROOT::VecOps;
-
-void rec_emcal_resolution_analysis(const char* input_fname = "rec_emcal_uniform_electrons.root")
-{
-  // 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();
-  ROOT::RDataFrame d0("events", input_fname);
-
-  // Number of Clusters
-  auto ncluster = [](const std::vector<eic::ClusterData>& evt) { return (int)evt.size(); };
-
-  // Cluster Energy [GeV]
-  auto clusterE = [](const std::vector<eic::ClusterData>& evt) {
-    std::vector<double> result;
-    for (const auto& i : evt)
-      result.push_back(i.energy);
-    return result;
-  };
-
-  // Thrown Energy [GeV]
-  auto E_thr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
-    std::vector<double> result;
-    result.push_back(TMath::Sqrt(input[2].psx * input[2].psx + input[2].psy * input[2].psy +
-                                 input[2].psz * input[2].psz + input[2].mass * input[2].mass));
-    return result;
-  };
-
-  // Energy Resolution = delta(E)/E
-  // delta(E) = (Erec - Ethr)
-  // E        = Ethr
-  auto E_res = [](const std::vector<double>& Erec, const std::vector<double>& Ethr) {
-    std::vector<double> result;
-    for (const auto& E1 : Ethr) {
-      for (const auto& E2 : Erec) {
-        result.push_back((E2 - E1) / E1);
-      }
-    }
-    return result;
-  };
-
-  // Define variables
-  auto d1 = d0.Define("ncluster", ncluster, {"EcalClusters"})
-                .Define("clusterE", clusterE, {"EcalClusters"})
-                .Define("E_thr", E_thr, {"mcparticles2"})
-                .Define("E_res", E_res, {"clusterE", "E_thr"});
-
-  // Define Cuts
-  // d2: Select Events with One Cluster
-  // d3: Select Events which determined their center of cluster not being on the edge of detector in addition to d2 cut
-  auto d2 = d1.Filter("ncluster==1");
-  auto d3 = d2.Filter(
-      [=](const std::vector<eic::ClusterData>& evt) {
-        for (const auto& i : evt) {
-          auto pos_x       = i.position.x;
-          auto pos_y       = i.position.y;
-          auto radial_dist = TMath::Sqrt(pos_x * pos_x + pos_y * pos_y);
-          if (radial_dist > 8.0 && radial_dist < 30.0)
-            return true;
-        }
-        return false;
-      },
-      {"EcalClusters"});
-
-  // Define histograms
-  auto hEresNC1 =
-      d2.Histo1D({"hEresNC1", "Energy Resolution CE vs Ethr w/ NC1; #DeltaE/E; Events", 100, -0.5, 0.5}, "E_res");
-  auto hEresNC1EGCUT = d3.Histo1D(
-      {"hEresNC1EGCUT", "Energy Resolution CE vs Ethr w/ NC1 EGCUT; #DeltaE/E; Events", 100, -0.5, 0.5}, "E_res");
-
-  // Energy Resolution Graphs
-  // Divided into serveral energy bin group
-  // Do gaussian fit to extract mean and sigma
-  int    nbin      = 12;                       // number of energy bin group
-  double maxEnergy = 30.0;                     // maximum of energy in dataset
-  double multiple  = maxEnergy / (double)nbin; // energy interval between energy bin group
-  double maximumE;                             // high end in each energy bin group
-  double minimumE;                             // low end in each energy bin group
-  double mean[nbin];                           // mean from gaussian fit
-  double sigma[nbin];                          // sigma from gaussian fit
-  double dmean[nbin];                          // error on mean
-  double dsigma[nbin];                         // error on sigma
-  double gEbin[nbin];                          // center of energy in each group
-  double invSqrtgEbin[nbin];                   // inverse sqrt center of energy in each group
-  double resolution[nbin];                     // resolution = sigma/mean
-  double resolutionError[nbin];                // error on resolution
-  double resolutionErrorSigma;                 // error on resolution sigma term
-  double resolutionErrorMean;                  // error on resolution mean term
-  double xerror[nbin];                         // set to 0.0 for all
-  // Loop over each energy bin group
-  // Define hitogram and do gaussian fit
-  // Save fit results; mean, sigma, error on mean & sigam
-  // Calculate resolution and resolution error
-  for (int n = 1; n <= nbin; n++) {
-    maximumE            = n * multiple;
-    minimumE            = maximumE - multiple;
-    gEbin[n - 1]        = (maximumE + minimumE) / 2.0;
-    invSqrtgEbin[n - 1] = 1.0 / TMath::Sqrt(gEbin[n - 1]);
-    // Define histogram
-    auto h = d3.Filter(
-                   [=](const std::vector<eic::ClusterData>& evt) {
-                     for (const auto& i : evt) {
-                       auto eng = i.energy;
-                       if (eng >= minimumE && eng < maximumE)
-                         return true;
-                     }
-                     return false;
-                   },
-                   {"EcalClusters"})
-                 .Histo1D({"h", "Energy Resolution; #DeltaE/E; Events", 100, -0.5, 0.5}, "E_res");
-    // Gaussian fit
-    h->Fit("gaus", "W","",h->GetMean() - 2.0 * h->GetRMS(),h->GetMean() + 2.0 * h->GetRMS());
-    // Save fit results
-    // how to access how many events in histogram: h->GetEntries();
-    if (h->GetFunction("gaus")) {
-      mean[n - 1]          = TMath::Abs(h->GetFunction("gaus")->GetParameter(1));
-      dmean[n - 1]         = h->GetFunction("gaus")->GetParError(1);
-      sigma[n - 1]         = h->GetFunction("gaus")->GetParameter(2);
-      dsigma[n - 1]        = h->GetFunction("gaus")->GetParError(2);
-      resolution[n - 1]    = sigma[n - 1] / mean[n - 1];
-      resolutionErrorSigma = dsigma[n - 1] / mean[n - 1];
-      resolutionErrorMean  = dmean[n - 1] * sigma[n - 1];
-      resolutionError[n - 1] =
-          TMath::Sqrt(resolutionErrorSigma * resolutionErrorSigma + resolutionErrorMean * resolutionErrorMean);
-      xerror[n - 1] = 0.0;
-    }
-    // Draw histogram, but
-    // Just save fit results
-    // TCanvas *c = new TCanvas("c", "c", 500, 500);
-    // h->GetYaxis()->SetTitleOffset(1.4);
-    // h->SetLineWidth(2);
-    // h->SetLineColor(kBlue);
-    // h->Fit("gaus","W,E");
-    // h->DrawClone();
-    // c->SaveAs("results/h.png");
-  }
-
-  // Generate graphs using gaussian fit results
-  // Cluster energy VS Simga
-  TGraphErrors* gCEVsSigmaNC1EGCUT = new TGraphErrors(nbin, gEbin, sigma, xerror, dsigma);
-  TCanvas*      c1                 = new TCanvas("c1", "c1", 500, 500);
-  gCEVsSigmaNC1EGCUT->SetTitle("Cluster Energy vs #sigma");
-  gCEVsSigmaNC1EGCUT->GetYaxis()->SetTitleOffset(1.4);
-  gCEVsSigmaNC1EGCUT->GetXaxis()->SetTitle("Cluster Energy [GeV]");
-  gCEVsSigmaNC1EGCUT->GetYaxis()->SetTitle("#sigma");
-  gCEVsSigmaNC1EGCUT->GetYaxis()->SetRangeUser(0.0, 0.05);
-  gCEVsSigmaNC1EGCUT->SetMarkerStyle(21);
-  gCEVsSigmaNC1EGCUT->DrawClone("AEP");
-  c1->SaveAs("results/emcal_electrons_gCEVsSigmaNC1EGCUT.png");
-  c1->SaveAs("results/emcal_electrons_gCEVsSigmaNC1EGCUT.pdf");
-  // Cluster energy VS Mean
-  TGraphErrors* gCEVsMeanNC1EGCUT = new TGraphErrors(nbin, gEbin, mean, xerror, dmean);
-  TCanvas*      c2                = new TCanvas("c2", "c2", 500, 500);
-  gCEVsMeanNC1EGCUT->SetTitle("Cluster Energy vs #mu");
-  gCEVsMeanNC1EGCUT->GetYaxis()->SetTitleOffset(1.4);
-  gCEVsMeanNC1EGCUT->GetXaxis()->SetTitle("Cluster Energy [GeV]");
-  gCEVsMeanNC1EGCUT->GetYaxis()->SetTitle("#mu");
-  gCEVsMeanNC1EGCUT->GetYaxis()->SetRangeUser(0.0, 0.1);
-  gCEVsMeanNC1EGCUT->SetMarkerStyle(21);
-  gCEVsMeanNC1EGCUT->DrawClone("AEP");
-  c2->SaveAs("results/emcal_electrons_gCEVsMeanNC1EGCUT.png");
-  c2->SaveAs("results/emcal_electrons_gCEVsMeanNC1EGCUT.pdf");
-  // Inverse of sqrt(cluster energy) VS resolution
-  TGraphErrors* gInvSqCEVsSigmaENC1EGCUT = new TGraphErrors(nbin, invSqrtgEbin, resolution, xerror, resolutionError);
-  TCanvas*      c3                       = new TCanvas("c3", "c3", 500, 500);
-  gInvSqCEVsSigmaENC1EGCUT->SetTitle("#sigma/E [%] vs 1/sqrt(Cluster Energy)");
-  gInvSqCEVsSigmaENC1EGCUT->GetYaxis()->SetTitleOffset(1.4);
-  gInvSqCEVsSigmaENC1EGCUT->GetXaxis()->SetTitle("1/sqrt(Cluster Energy)");
-  gInvSqCEVsSigmaENC1EGCUT->GetYaxis()->SetTitle("#sigma/E [%]");
-  gInvSqCEVsSigmaENC1EGCUT->GetYaxis()->SetRangeUser(0.0, 1.5);
-  gInvSqCEVsSigmaENC1EGCUT->GetXaxis()->SetRangeUser(0.0, 1.5);
-  gInvSqCEVsSigmaENC1EGCUT->SetMarkerStyle(21);
-  gInvSqCEVsSigmaENC1EGCUT->DrawClone("AEP");
-  c3->SaveAs("results/emcal_electrons_gInvSqCEVsSimgaENC1EGCUT.png");
-  c3->SaveAs("results/emcal_electrons_gInvSqCEVsSigmaENC1EGCUT.pdf");
-
-  // Fit Function
-  // sigma/E = S/sqrt(E) + N/E + C; added in quadrature
-  // S: Stochastic term
-  // N: Niose term
-  // C: Constant term
-  TF1* f1 = new TF1("f1", "sqrt([0]*[0] + pow([1]/sqrt(x),2) + pow([2]/x,2))", 0.5, 30.0);
-  f1->SetParNames("C", "S", "N");      // param names
-  f1->SetParameters(0.30, 2.80, 0.12); // initial param values
-  f1->SetLineWidth(2);
-  f1->SetLineColor(kRed);
-
-  // Cluster energy VS resolution
-  TGraphErrors* gCEVsSigmaENC1EGCUT = new TGraphErrors(nbin, gEbin, resolution, xerror, resolutionError);
-  TCanvas*      c4                  = new TCanvas("c4", "c4", 500, 500);
-  gCEVsSigmaENC1EGCUT->SetTitle("Cluster Energy vs #sigma/E [%]");
-  gCEVsSigmaENC1EGCUT->GetYaxis()->SetTitleOffset(1.4);
-  gCEVsSigmaENC1EGCUT->GetXaxis()->SetTitle("Cluster Energy [GeV]");
-  gCEVsSigmaENC1EGCUT->GetYaxis()->SetTitle("#sigma/E [%]");
-  gCEVsSigmaENC1EGCUT->GetYaxis()->SetRangeUser(0.0, 1.5);
-  gCEVsSigmaENC1EGCUT->SetMarkerStyle(21);
-  gCEVsSigmaENC1EGCUT->Fit("f1", "W");
-  gCEVsSigmaENC1EGCUT->DrawClone("AEP");
-  c4->SaveAs("results/emcal_electrons_gCEVsSimgaENC1EGCUT.png");
-  c4->SaveAs("results/emcal_electrons_gCEVsSigmaENC1EGCUT.pdf");
-
-  // Event Counts
-  auto nevents_thrown      = d1.Count();
-  auto nevents_cluster1    = d2.Count();
-  auto nevents_cluster1cut = d3.Count();
-
-  std::cout << "Number of Events: " << (*nevents_cluster1) << " / " << (*nevents_thrown)
-            << " = single cluster events/all \n";
-  std::cout << "Number of Events: " << (*nevents_cluster1cut) << " / " << (*nevents_thrown)
-            << " = single cluster events that passed edge cut/all \n";
-}
diff --git a/benchmarks/ecal/scripts/draw_clusters.py b/benchmarks/ecal/scripts/draw_clusters.py
index f4c81896..74c7fe22 100644
--- a/benchmarks/ecal/scripts/draw_clusters.py
+++ b/benchmarks/ecal/scripts/draw_clusters.py
@@ -26,7 +26,6 @@ def load_root_macros(arg_macros):
 def flatten_collection(rdf, collection, cols=None):
     if not cols:
         cols  = [str(c) for c in rdf.GetColumnNames() if str(c).startswith('{}.'.format(collection))]
-        cols += [str(c) for c in rdf.GetColumnNames() if str(c).startswith('{}Info.'.format(collection))]
     else:
         cols = ['{}.{}'.format(collection, c) for c in cols]
     if not cols:
@@ -93,8 +92,6 @@ if __name__ == '__main__':
     ]
     for coll in [c.strip() for c in args.coll.split(',')]:
         df = flatten_collection(rdf_rec, coll)
-        # Get eta from info table
-        df[coll + '.eta'] = df[coll + 'Info.eta']
         fig, axs = plt.subplots(2, 2, figsize=(12, 8), dpi=160)
         ncl = df.groupby('event')['{}.ID.value'.format(coll)].nunique().values
         axs[0][0].hist(ncl, weights=np.repeat(1./float(ncl.shape[0]), ncl.shape[0]),
diff --git a/benchmarks/imaging_ecal/options/hybrid_cluster.py b/benchmarks/imaging_ecal/options/hybrid_cluster.py
index 21fef284..295be391 100644
--- a/benchmarks/imaging_ecal/options/hybrid_cluster.py
+++ b/benchmarks/imaging_ecal/options/hybrid_cluster.py
@@ -91,7 +91,7 @@ clusterreco = ImagingClusterReco('imcal_clreco',
         inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection,
         outputLayerCollection='EcalBarrelImagingClustersLayers',
         outputClusterCollection='EcalBarrelImagingClusters',
-        outputInfoCollection='EcalBarrelImagingClustersInfo',
+        mcHits="EcalBarrelHits",
         samplingFraction=kwargs['img_sf'])
 
 # scfi layers
@@ -138,7 +138,7 @@ scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
        inputHitCollection=scfi_barrel_cl.inputHitCollection,
        inputProtoClusterCollection=scfi_barrel_cl.outputProtoClusterCollection,
        outputClusterCollection="EcalBarrelScFiClusters",
-       outputInfoCollection="EcalBarrelScFiClustersInfo",
+        mcHits="EcalBarrelScFiHits",
        logWeightBase=6.2,
        samplingFraction=kwargs['scfi_sf'])
 
diff --git a/benchmarks/imaging_ecal/options/imaging_2dcluster.py b/benchmarks/imaging_ecal/options/imaging_2dcluster.py
index cf0a1794..9374a8f7 100644
--- a/benchmarks/imaging_ecal/options/imaging_2dcluster.py
+++ b/benchmarks/imaging_ecal/options/imaging_2dcluster.py
@@ -85,7 +85,7 @@ imcal_barrel_clreco = RecoCoG('imcal_barrel_clreco',
         inputHitCollection=imcal_barrel_cl.inputHitCollection,
         inputProtoClusterCollection = imcal_barrel_cl.outputProtoClusterCollection,
         outputClusterCollection='EcalBarrelImagingClusters',
-        outputInfoCollection='EcalBarrelImagingClustersInfo',
+        mcHits="EcalBarrelHits",
         logWeightBase=6.2,
         samplingFraction=kwargs['img_sf'])
 
diff --git a/benchmarks/imaging_ecal/options/imaging_topocluster.py b/benchmarks/imaging_ecal/options/imaging_topocluster.py
index 2ce1fc9a..a42c3a10 100644
--- a/benchmarks/imaging_ecal/options/imaging_topocluster.py
+++ b/benchmarks/imaging_ecal/options/imaging_topocluster.py
@@ -85,7 +85,7 @@ clusterreco = ImagingClusterReco("imcal_clreco",
         inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection,
         outputLayerCollection="EcalBarrelImagingClustersLayers",
         outputClusterCollection="EcalBarrelImagingClusters",
-        outputInfoCollection="EcalBarrelImagingClustersInfo",
+        mcHits="EcalBarrelHits",
         samplingFraction=sf)
 
 out.outputCommands = ["keep *"]
diff --git a/benchmarks/imaging_ecal/options/scfi_cluster.py b/benchmarks/imaging_ecal/options/scfi_cluster.py
index 1c0def2b..741d124a 100644
--- a/benchmarks/imaging_ecal/options/scfi_cluster.py
+++ b/benchmarks/imaging_ecal/options/scfi_cluster.py
@@ -82,7 +82,7 @@ scfi_barrel_clreco = RecoCoG("scfi_barrel_clreco",
        inputHitCollection=scfi_barrel_cl.inputHitCollection,
        inputProtoClusterCollection=scfi_barrel.outputProtoClusterCollection,
        outputClusterCollection="EcalBarrelScFiClusters",
-       outputInfoCollection="EcalBarrelScFiClustersInfo",
+       mcHits="EcalBarrelScFiHits",
        logWeightBase=6.2,
        samplingFraction=kwargs['sf'])
 
diff --git a/benchmarks/imaging_ecal/scripts/aggregator_occupancy.py b/benchmarks/imaging_ecal/scripts/aggregator_occupancy.py
index d10661b9..62406355 100644
--- a/benchmarks/imaging_ecal/scripts/aggregator_occupancy.py
+++ b/benchmarks/imaging_ecal/scripts/aggregator_occupancy.py
@@ -29,7 +29,6 @@ def load_root_macros(arg_macros):
 def flatten_collection(rdf, collection, cols=None):
     if not cols:
         cols = [str(c) for c in rdf.GetColumnNames() if str(c).startswith('{}.'.format(collection))]
-        cols += [str(c) for c in rdf.GetColumnNames() if str(c).startswith('{}Info.'.format(collection))]
     else:
         cols = ['{}.{}'.format(collection, c) for c in cols]
     if not cols:
-- 
GitLab