diff --git a/benchmarks/clustering/barrel_clusters.sh b/benchmarks/clustering/barrel_clusters.sh
index 4d671c8b225104603b6d88cd00040940571d5f58..a34bd51b57b45117e2b2e85114f9f3c34f1660fc 100644
--- a/benchmarks/clustering/barrel_clusters.sh
+++ b/benchmarks/clustering/barrel_clusters.sh
@@ -19,7 +19,7 @@ if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
fi
if [[ ! -n "${CB_EMCAL_SAMP_FRAC}" ]] ; then
- export CB_EMCAL_SAMP_FRAC=0.014
+ export CB_EMCAL_SAMP_FRAC=1.0
fi
export JUGGLER_FILE_NAME_TAG="barrel_clusters"
@@ -47,7 +47,7 @@ if [[ "$?" -ne "0" ]] ; then
fi
# Need to figure out how to pass file name to juggler from the commandline
-gaudirun.py benchmarks/clustering/options/fullcalo_clustering.py
+gaudirun.py benchmarks/clustering/options/hcal_clustering.py
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running juggler"
exit 1
diff --git a/benchmarks/clustering/options/hcal_clustering.py b/benchmarks/clustering/options/hcal_clustering.py
new file mode 100644
index 0000000000000000000000000000000000000000..bdd8a2af679be9cb192ff4103f0e84714b839cd6
--- /dev/null
+++ b/benchmarks/clustering/options/hcal_clustering.py
@@ -0,0 +1,136 @@
+from Gaudi.Configuration import *
+import os
+
+from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
+from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+from GaudiKernel import SystemOfUnits as units
+
+detector_name = "topside"
+if "JUGGLER_DETECTOR" in os.environ :
+ detector_name = str(os.environ["JUGGLER_DETECTOR"])
+
+if "JUGGLER_DETECTOR_PATH" in os.environ :
+ detector_name = str(os.environ["JUGGLER_DETECTOR_PATH"]) + "/" + detector_name
+
+# get sampling fraction from system environment variable, 1.0 by default
+sf = 1.0
+if "CB_EMCAL_SAMP_FRAC" in os.environ :
+ sf = str(os.environ["CB_EMCAL_SAMP_FRAC"])
+
+# todo add checks
+input_sim_file = str(os.environ["JUGGLER_SIM_FILE"])
+output_rec_file = str(os.environ["JUGGLER_REC_FILE"])
+n_events = str(os.environ["JUGGLER_N_EVENTS"])
+
+geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(detector_name)])
+podioevent = EICDataSvc("EventDataSvc", inputs=[input_sim_file])
+
+from Configurables import PodioInput
+from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
+from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
+
+# from Configurables import Jug__Digi__HadronicCalDigi as HadCalorimeterDigi
+from Configurables import Jug__Digi__CalorimeterHitDigi as HadCalorimeterDigi
+from Configurables import Jug__Digi__EMCalorimeterDigi as EMCalorimeterDigi
+from Configurables import Jug__Digi__CrystalEndcapsDigi as CrystalEndcapsDigi
+
+from Configurables import Jug__Reco__CrystalEndcapsReco as CrystalEndcapsReco
+from Configurables import Jug__Reco__EMCalReconstruction as EMCalReconstruction
+# from Configurables import Jug__Reco__HCalReconstruction as HCalReconstruction
+from Configurables import Jug__Reco__EcalTungstenSamplingReco as HCalReconstruction
+
+# from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
+from Configurables import Jug__Reco__TopologicalCellCluster as TopoCluster
+from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
+from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
+
+podioinput = PodioInput("PodioReader", collections=["mcparticles","HcalBarrelHits"], OutputLevel=DEBUG)
+
+## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
+copier = MCCopier("MCCopier",
+ inputCollection="mcparticles",
+ outputCollection="mcparticles2")
+
+##raw hits
+
+
+calcopier = CalCopier("CalCopier",
+ inputCollection="HcalBarrelHits",
+ outputCollection="HcalBarrelHits2")
+
+#digitized hits
+
+hcaldigi = HadCalorimeterDigi("hcal_barrel_digi",
+ inputHitCollection="HcalBarrelHits",
+ outputHitCollection="RawHcalBarrelHits",
+ inputEnergyUnit=units.GeV,
+ inputTimeUnit=units.ns,
+ energyResolutions=[0.07, 0., 0.],
+ dynamicRangeADC=2.*units.GeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=10,
+ OutputLevel=DEBUG
+ )
+
+#reconstructed hits
+
+hcal_reco = HCalReconstruction("hcal_reco",
+ inputHitCollection="RawHcalBarrelHits",
+ outputHitCollection="RecHcalBarrelHits",
+ dynamicRangeADC=2.*units.GeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=10,
+ thresholdFactor=5.0,
+ OutputLevel=DEBUG)
+
+#clusters
+
+#hcal_barrel_cluster = TopoCluster("hcal_barrel_cluster",
+# inputHitCollection="RecHcalBarrelHits",
+# outputClusterCollection="HcalBarrelClusters",
+# adjLayerDiff=2,
+# localRanges=[20.*units.mm, 20.*units.mm], # local x, y for hits at the same layer
+# adjLayerRanges=[0.02, 0.02], # eta, phi for different layers, the same sector
+# adjSectorDist=5.*units.cm, # different sector
+# minClusterCenterEdep=10.*units.MeV,
+# readoutClass="HcalBarrelHits", # readout class name to get layer/sector id
+# layerField="layer",
+# sectorField="module",
+# OutputLevel=DEBUG
+# )
+
+
+hcal_barrel_cluster = IslandCluster("hcal_barrel_cluster",
+ inputHitCollection="RecHcalBarrelHits",
+ outputClusterCollection="HcalBarrelClusters",
+ minClusterCenterEdep=30*units.MeV,
+ groupRange=2.0,
+ OutputLevel=DEBUG)
+
+
+# finalizing clustering (center-of-gravity step)
+
+hcal_barrel_clusterreco = RecoCoG("hcal_barrel_clusterreco",
+ clusterCollection="HcalBarrelClusters",
+ logWeightBase=6.2,
+ samplingFraction=sf,
+ OutputLevel=DEBUG)
+
+out = PodioOutput("out", filename=output_rec_file)
+
+out.outputCommands = ["keep *"]
+
+ApplicationMgr(
+ TopAlg = [podioinput, copier, calcopier,
+ hcaldigi,
+ hcal_reco,
+ hcal_barrel_cluster,
+ hcal_barrel_clusterreco,
+ out],
+ EvtSel = 'NONE',
+ EvtMax = n_events,
+ ExtSvc = [podioevent],
+ OutputLevel=DEBUG
+ )
diff --git a/benchmarks/clustering/scripts/gen_central_pions.cxx b/benchmarks/clustering/scripts/gen_central_pions.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..80aa4fcba96045ced2cdb707d1c0ad4447da332f
--- /dev/null
+++ b/benchmarks/clustering/scripts/gen_central_pions.cxx
@@ -0,0 +1,80 @@
+#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;
+
+/** Generate electrons in the central region.
+ * This is for testing detectors in the "barrel" region.
+ */
+void gen_central_pions(int n_events = 100,
+ const char* out_fname = "central_pions.hepmc")
+{
+ double cos_theta_min = std::cos( 60.0*(M_PI/180.0));
+ double cos_theta_max = std::cos(120.0*(M_PI/180.0));
+
+ WriterAscii hepmc_output(out_fname);
+ int events_parsed = 0;
+ GenEvent evt(Units::GEV, Units::MM);
+
+ // Random number generator
+ TRandom *r1 = new TRandom();
+
+ 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(100.0, 100.1);
+ Double_t phi = r1->Uniform(0.0, 0.25 * M_PI);
+ Double_t costh = r1->Uniform(cos_theta_min, cos_theta_max);
+ Double_t th = std::acos(costh);
+ Double_t px = p * std::cos(phi) * std::sin(th);
+ Double_t py = p * std::sin(phi) * std::sin(th);
+ Double_t pz = p * std::cos(th);
+ // Generates random vectors, uniformly distributed over the surface of a
+ // sphere of given radius, in this case momentum.
+ // r1->Sphere(px, py, pz, p);
+
+ std::cout << std::sqrt(px*px + py*py + pz*pz) - p << " is zero? \n";
+
+ // 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.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;
+}