diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index a031ed1fc421ed420f91589d922661e61b2b7bb2..beb4aa328fff9504acedf655e1647b810ad0d118 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -49,7 +49,7 @@ include:
final_report:
stage: finish
- needs: ["ecal_1_emcal_electrons", "tracking_central_electrons"]
+ needs: ["ecal_1_emcal_crystal_electrons", "tracking_central_electrons"]
script:
- mkdir -p results/views && cd results/views && bash ../../bin/download_views
- echo "It was a success!"
diff --git a/benchmarks/clustering/options/calorimeter_clustering.py b/benchmarks/clustering/options/calorimeter_clustering.py
index a45e642a7cc4e694a14dce0cc6cf7aba21a27f3f..bb98b16b24476ed8beee0a822459ff7023ef8e84 100644
--- a/benchmarks/clustering/options/calorimeter_clustering.py
+++ b/benchmarks/clustering/options/calorimeter_clustering.py
@@ -35,7 +35,7 @@ from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
-podioinput = PodioInput("PodioReader", collections=["mcparticles","CrystalEcalHits","EcalBarrelHits"], OutputLevel=DEBUG)
+podioinput = PodioInput("PodioReader", collections=["mcparticles","CrystalEcalHits","EcalBarrelAstroPixHits"], OutputLevel=DEBUG)
## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
copier = MCCopier("MCCopier",
@@ -49,8 +49,8 @@ emcaldigi = CrystalEndcapsDigi("ecal_digi",
inputHitCollection="CrystalEcalHits",
outputHitCollection="RawDigiEcalHits")
ecdigi = EMCalorimeterDigi("ec_barrel_digi",
- inputHitCollection="EcalBarrelHits",
- outputHitCollection="RawEcalBarrelHits")
+ inputHitCollection="EcalBarrelAstroPixHits",
+ outputHitCollection="RawEcalBarrelAstroPixHits")
crystal_ec_reco = CrystalEndcapsReco("crystal_ec_reco",
inputHitCollection="RawDigiEcalHits",
@@ -58,15 +58,15 @@ crystal_ec_reco = CrystalEndcapsReco("crystal_ec_reco",
minModuleEdep=1.0*units.MeV)
ecal_reco = EMCalReconstruction("ecal_reco",
- inputHitCollection="RawEcalBarrelHits",
- outputHitCollection="RecEcalBarrelHits",
+ inputHitCollection="RawEcalBarrelAstroPixHits",
+ outputHitCollection="RecEcalBarrelAstroPixHits",
samplingFraction=0.25,
minModuleEdep=0.0*units.MeV)
ec_barrel_cluster = IslandCluster("ec_barrel_cluster",
- inputHitCollection="RecEcalBarrelHits",
- outputClusterCollection="EcalBarrelClusters",
- splitHitCollection="splitEcalBarrelHitCollection",
+ inputHitCollection="RecEcalBarrelAstroPixHits",
+ outputClusterCollection="EcalBarrelAstroPixClusters",
+ splitHitCollection="splitEcalBarrelAstroPixHitCollection",
minClusterCenterEdep=1*units.MeV,
groupRange=2.0,
OutputLevel=DEBUG)
@@ -78,7 +78,7 @@ crystal_ec_cluster = IslandCluster("crystal_ec_cluster",
OutputLevel=DEBUG)
simple_cluster = SimpleClustering("simple_cluster",
- inputHitCollection="RecEcalBarrelHits",
+ inputHitCollection="RecEcalBarrelAstroPixHits",
outputClusters="SimpleClusters",
OutputLevel=DEBUG)
diff --git a/benchmarks/ecal/config.yml b/benchmarks/ecal/config.yml
index ffdfdd73fc39948ff26e9c64e4af3a3acbf53d61..814f77c70b99a68203626a8144c7e3a29d6df0c3 100644
--- a/benchmarks/ecal/config.yml
+++ b/benchmarks/ecal/config.yml
@@ -1,4 +1,4 @@
-ecal_1_emcal_electrons:
+ecal_1_emcal_crystal_electrons:
image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG
stage: process
needs: ["detector"]
@@ -10,7 +10,7 @@ ecal_1_emcal_electrons:
script:
- bash benchmarks/ecal/emcal_electrons.sh
-ecal_1_emcal_pi0s:
+ecal_2_emcal_crystal_pi0s:
image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG
needs: ["detector"]
timeout: 12 hours
@@ -22,7 +22,32 @@ ecal_1_emcal_pi0s:
script:
- bash benchmarks/ecal/emcal_pi0s.sh
-emcal_barrel_electrons:
+ecal_3_emcal_barrel_electrons:
+ image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG
+ needs: ["detector"]
+ timeout: 48 hours
+ artifacts:
+ expire_in: 20 weeks
+ paths:
+ - results/
+ stage: run
+ script:
+ - bash benchmarks/ecal/run_emcal_barrel_electrons.sh
+
+ecal_4_emcal_barrel_pions:
+ image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG
+ needs: ["detector"]
+ timeout: 48 hours
+ artifacts:
+ expire_in: 20 weeks
+ paths:
+ - results/
+ stage: run
+ script:
+ - bash benchmarks/ecal/run_emcal_barrel_pions.sh
+
+## Not fully yet...
+full_emcal_barrel_electrons:
image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG
needs: ["detector"]
timeout: 48 hours
diff --git a/benchmarks/ecal/options/emcal_barrel_reco.py b/benchmarks/ecal/options/emcal_barrel_reco.py
new file mode 100644
index 0000000000000000000000000000000000000000..f43b1d3e7c57b97bb08319d90dd4401c4b2e17ce
--- /dev/null
+++ b/benchmarks/ecal/options/emcal_barrel_reco.py
@@ -0,0 +1,115 @@
+#######################################
+# EMCAL Barrel detector Reconstruction
+# J.KIM 04/02/2021
+#######################################
+
+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
+
+input_sim_file = "jug_input.root"
+if "JUGGLER_SIM_FILE" in os.environ :
+ input_sim_file = str(os.environ["JUGGLER_SIM_FILE"])
+else :
+ print(" ERROR : JUGGLER_SIM_FILE not set" )
+
+output_rec_file = "jug_rec.root"
+if "JUGGLER_REC_FILE" in os.environ :
+ output_rec_file = str(os.environ["JUGGLER_REC_FILE"])
+else :
+ print(" ERROR : JUGGLER_REC_FILE not set" )
+
+n_events = 100
+if "JUGGLER_N_EVENTS" in os.environ :
+ n_events = str(os.environ["JUGGLER_N_EVENTS"])
+
+geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(detector_name)])
+podioevent = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=DEBUG)
+
+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__EcalTungstenSamplingDigi as EcalTungstenSamplingDigi
+
+from Configurables import Jug__Reco__EcalTungstenSamplingReco as EcalTungstenSamplingReco
+from Configurables import Jug__Reco__SamplingECalHitsMerger as SamplingECalHitsMerger
+from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
+from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
+
+podioinput = PodioInput("PodioReader", collections=["mcparticles","EcalBarrelAstroPixHits"], OutputLevel=DEBUG)
+
+# Thrown Information
+copier = MCCopier("MCCopier",
+ inputCollection="mcparticles",
+ outputCollection="mcparticles2",
+ OutputLevel=DEBUG)
+# Geant4 Information
+embarrelcopier = CalCopier("CalBarrelCopier",
+ inputCollection="EcalBarrelAstroPixHits",
+ outputCollection="EcalBarrelAstroPixHits2",
+ OutputLevel=DEBUG)
+# Digitization
+embarreldigi = EcalTungstenSamplingDigi("ecal_barrel_digi",
+ inputHitCollection="EcalBarrelAstroPixHits",
+ outputHitCollection="RawEcalBarrelAstroPixHits",
+ inputEnergyUnit=units.GeV,
+ inputTimeUnit=units.ns,
+ dynamicRangeADC=700*units.keV,
+ energyResolutions=[0., 0.02, 0.],
+ pedestalSigma=40,
+ OutputLevel=DEBUG)
+# Reconstruction
+embarrelreco = EcalTungstenSamplingReco("ecal_barrel_reco",
+ inputHitCollection="RawEcalBarrelAstroPixHits",
+ outputHitCollection="RecoEcalBarrelAstroPixHits",
+ dynamicRangeADC=700*units.keV,
+ pedestalSigma=40,
+ OutputLevel=DEBUG)
+# 2D+1 Clusterings
+# readout id definition for barrel ecal
+# <id>system:8,barrel:3,module:4,layer:10,slice:5,x:32:-16,y:-16</id>
+# xy_merger sum layers/slices, masking (8+3+4, 8+3+4+5+10-1)
+embarrelxymerger = SamplingECalHitsMerger("ecal_barrel_xy_merger",
+ cellIDMaskRanges=[(15, 29)],
+ inputHitCollection="RecoEcalBarrelAstroPixHits",
+ outputHitCollection="RecoEcalBarrelAstroPixHitsXY")
+# xy_merger sum modules, masking (8+3+4+5+10, 8+3+4+5+10+32-1)
+embarrelzmerger = SamplingECalHitsMerger("ecal_barrel_z_merger",
+ cellIDMaskRanges=[(30, 61)],
+ inputHitCollection="RecoEcalBarrelAstroPixHits",
+ outputHitCollection="RecoEcalBarrelAstroPixHitsZ")
+# Clustering
+embarrelcluster = IslandCluster("ecal_barrel_cluster",
+ inputHitCollection="RecoEcalBarrelAstroPixHitsXY",
+ outputClusterCollection="EcalBarrelAstroPixClusters",
+ minClusterCenterEdep=0.5*units.MeV,
+ splitCluster=False,
+ groupRanges=[2.0*units.cm, 2.0*units.cm, 2.0*units.cm])
+# Reconstruct the cluster with Center of Gravity method
+embarrelclusterreco = RecoCoG("ecal_barrel_clusterreco",
+ clusterCollection="EcalBarrelAstroPixClusters",
+ logWeightBase=6.2)
+
+out = PodioOutput("out", filename=output_rec_file)
+
+out.outputCommands = ["keep *"]
+
+ApplicationMgr(
+ TopAlg = [podioinput, copier, embarrelcopier, embarreldigi,
+ embarrelreco, embarrelxymerger, embarrelzmerger, embarrelcluster, embarrelclusterreco, out],
+ EvtSel = 'NONE',
+ EvtMax = n_events,
+ ExtSvc = [podioevent],
+ OutputLevel=DEBUG
+ )
diff --git a/benchmarks/ecal/options/full_em_calorimeter_reco.py b/benchmarks/ecal/options/full_em_calorimeter_reco.py
index a61364b0234ac3c3958351e360829cb77e0ee391..a0c6d0f2e1633304db3fbf45826f31e57ec76119 100644
--- a/benchmarks/ecal/options/full_em_calorimeter_reco.py
+++ b/benchmarks/ecal/options/full_em_calorimeter_reco.py
@@ -42,7 +42,7 @@ from Configurables import Jug__Reco__CrystalEndcapsReco as CrystalEndcapsReco
from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
-podioinput = PodioInput("PodioReader", collections=["mcparticles","CrystalEcalHits","EcalBarrelHits","EcalEndcapHits"], OutputLevel=DEBUG)
+podioinput = PodioInput("PodioReader", collections=["mcparticles","CrystalEcalHits","EcalBarrelAstroPixHits","EcalEndcapHits"], OutputLevel=DEBUG)
## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
copier = MCCopier("MCCopier",
@@ -55,10 +55,19 @@ calcopier = CalCopier("CalCopier",
outputCollection="CrystalEcalHits2",
OutputLevel=DEBUG)
-embarreldigi = EcalTungstenSamplingDigi("ec_barrel_digi",
- inputHitCollection="EcalBarrelHits",
- outputHitCollection="RawEcalBarrelHits",
- energyResolution=0.11,
+embarrelcopier = CalCopier("CalBarrelCopier",
+ inputCollection="EcalBarrelAstroPixHits",
+ outputCollection="EcalBarrelAstroPixHits2",
+ OutputLevel=DEBUG)
+
+embarreldigi = EcalTungstenSamplingDigi("ecal_barrel_digi",
+ inputHitCollection="EcalBarrelAstroPixHits",
+ outputHitCollection="RawEcalBarrelAstroPixHits",
+ inputEnergyUnit=units.GeV,
+ inputTimeUnit=units.ns,
+ dynamicRangeADC=700*units.keV,
+ energyResolutions=[0., 0.02, 0.],
+ pedestalSigma=40,
OutputLevel=DEBUG)
emendcapdigi = EcalTungstenSamplingDigi("ec_endcap_digi",
@@ -92,7 +101,7 @@ out = PodioOutput("out", filename=output_rec_file)
out.outputCommands = ["keep *"]
ApplicationMgr(
- TopAlg = [podioinput, copier, calcopier,
+ TopAlg = [podioinput, copier, calcopier, embarrelcopier,
embarreldigi, emendcapdigi, emcaldigi, emcalreco, emcalcluster, clusterreco, out],
EvtSel = 'NONE',
EvtMax = n_events,
diff --git a/benchmarks/ecal/run_emcal_barrel_electrons.sh b/benchmarks/ecal/run_emcal_barrel_electrons.sh
new file mode 100755
index 0000000000000000000000000000000000000000..9b6edf1fb0552c346c21209369e424fc0205fb6d
--- /dev/null
+++ b/benchmarks/ecal/run_emcal_barrel_electrons.sh
@@ -0,0 +1,89 @@
+#!/bin/bash
+
+./util/print_env.sh
+
+## To run the reconstruction, we need the following global variables:
+## - JUGGLER_INSTALL_PREFIX: Install prefix for Juggler (simu/recon)
+## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark
+## - JUGGLER_DETECTOR_VERSION: the detector package we want to use for this benchmark
+## - DETECTOR_PATH: full path to the detector definitions
+##
+## You can ready options/env.sh for more in-depth explanations of the variables
+## and how they can be controlled.
+source options/env.sh
+export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH}
+
+if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
+ export JUGGLER_N_EVENTS=1000
+fi
+
+if [[ ! -n "${E_start}" ]] ; then
+ export E_start=5.0
+fi
+
+if [[ ! -n "${E_end}" ]] ; then
+ export E_end=5.0
+fi
+
+export JUGGLER_FILE_NAME_TAG="emcal_barrel_uniform_electrons"
+export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
+
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
+
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
+
+# Generate the input events
+root -b -q "benchmarks/ecal/scripts/emcal_barrel_electrons.cxx(${JUGGLER_N_EVENTS}, ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running script: generating input events"
+ exit 1
+fi
+# Plot the input events
+root -b -q "benchmarks/ecal/scripts/emcal_barrel_electrons_reader.cxx(${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running script: plotting input events"
+ exit 1
+fi
+
+# Run geant4 simulations
+npsim --runType batch \
+ -v WARNING \
+ --part.minimalKineticEnergy 0.5*GeV \
+ --numberOfEvents ${JUGGLER_N_EVENTS} \
+ --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
+ --inputFiles ${JUGGLER_FILE_NAME_TAG}.hepmc \
+ --outputFile ${JUGGLER_SIM_FILE}
+
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running npdet"
+ exit 1
+fi
+
+# Run Juggler
+xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+ gaudirun.py benchmarks/ecal/options/emcal_barrel_reco.py
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running juggler"
+ exit 1
+fi
+
+# Directory for plots
+mkdir -p results
+
+# Run analysis script
+root -b -q "benchmarks/ecal/scripts/emcal_barrel_electrons_analysis.cxx(\"${JUGGLER_REC_FILE}\")"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running analysis script"
+ exit 1
+fi
+
+root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}")
+if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
+ # file must be less than 10 MB to upload
+ if [[ "${root_filesize}" -lt "10000000" ]] ; then
+ cp ${JUGGLER_REC_FILE} results/.
+ fi
+fi
+
diff --git a/benchmarks/ecal/run_emcal_barrel_pions.sh b/benchmarks/ecal/run_emcal_barrel_pions.sh
new file mode 100755
index 0000000000000000000000000000000000000000..5115a63d9cd097b7de6819db6453ac3f97aea8b8
--- /dev/null
+++ b/benchmarks/ecal/run_emcal_barrel_pions.sh
@@ -0,0 +1,89 @@
+#!/bin/bash
+
+./util/print_env.sh
+
+## To run the reconstruction, we need the following global variables:
+## - JUGGLER_INSTALL_PREFIX: Install prefix for Juggler (simu/recon)
+## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark
+## - JUGGLER_DETECTOR_VERSION: the detector package we want to use for this benchmark
+## - DETECTOR_PATH: full path to the detector definitions
+##
+## You can ready options/env.sh for more in-depth explanations of the variables
+## and how they can be controlled.
+source options/env.sh
+export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH}
+
+if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
+ export JUGGLER_N_EVENTS=1000
+fi
+
+if [[ ! -n "${E_start}" ]] ; then
+ export E_start=5.0
+fi
+
+if [[ ! -n "${E_end}" ]] ; then
+ export E_end=5.0
+fi
+
+export JUGGLER_FILE_NAME_TAG="emcal_barrel_uniform_pions"
+export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
+
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
+
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
+
+# Generate the input events
+root -b -q "benchmarks/ecal/scripts/emcal_barrel_pions.cxx(${JUGGLER_N_EVENTS}, ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running script: generating input events"
+ exit 1
+fi
+# Plot the input events
+root -b -q "benchmarks/ecal/scripts/emcal_barrel_pions_reader.cxx(${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running script: plotting input events"
+ exit 1
+fi
+
+# Run geant4 simulations
+npsim --runType batch \
+ -v WARNING \
+ --part.minimalKineticEnergy 0.5*GeV \
+ --numberOfEvents ${JUGGLER_N_EVENTS} \
+ --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
+ --inputFiles ${JUGGLER_FILE_NAME_TAG}.hepmc \
+ --outputFile ${JUGGLER_SIM_FILE}
+
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running npdet"
+ exit 1
+fi
+
+# Run Juggler
+xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+ gaudirun.py benchmarks/ecal/options/emcal_barrel_reco.py
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running juggler"
+ exit 1
+fi
+
+# Directory for plots
+mkdir -p results
+
+# Run analysis script
+root -b -q "benchmarks/ecal/scripts/emcal_barrel_pions_analysis.cxx(\"${JUGGLER_REC_FILE}\")"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running analysis script"
+ exit 1
+fi
+
+root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}")
+if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
+ # file must be less than 10 MB to upload
+ if [[ "${root_filesize}" -lt "10000000" ]] ; then
+ cp ${JUGGLER_REC_FILE} results/.
+ fi
+fi
+
diff --git a/benchmarks/ecal/scripts/emcal_barrel_electrons.cxx b/benchmarks/ecal/scripts/emcal_barrel_electrons.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..12a78da363d28a2725f022685a8306aa8532ec2f
--- /dev/null
+++ b/benchmarks/ecal/scripts/emcal_barrel_electrons.cxx
@@ -0,0 +1,80 @@
+//////////////////////////////////////////////////////////////
+// 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/emcal_barrel_electrons_analysis.cxx b/benchmarks/ecal/scripts/emcal_barrel_electrons_analysis.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..f516ec63bc122d179cd152d2b3285291b04cf7dc
--- /dev/null
+++ b/benchmarks/ecal/scripts/emcal_barrel_electrons_analysis.cxx
@@ -0,0 +1,170 @@
+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/CalorimeterHitCollection.h"
+#include "eicd/RawCalorimeterHitCollection.h"
+#include "eicd/RawCalorimeterHitData.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 "TF1.h"
+#include "TH1D.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, {"RecoEcalBarrelAstroPixHitsXY"})
+ .Define("ErecZ", ErecZ, {"RecoEcalBarrelAstroPixHitsZ"})
+ .Define("ncluster", ncluster, {"EcalBarrelAstroPixClusters"})
+ .Define("Ecluster", Ecluster, {"EcalBarrelAstroPixClusters"})
+ .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/emal_electrons_ErecZ.png");
+ c3->SaveAs("results/emal_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);
+ 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/emcal_barrel_electrons_reader.cxx b/benchmarks/ecal/scripts/emcal_barrel_electrons_reader.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..75c2dae9c9218f036bbca7dd65f65e85bd67eca2
--- /dev/null
+++ b/benchmarks/ecal/scripts/emcal_barrel_electrons_reader.cxx
@@ -0,0 +1,125 @@
+//////////////////////////
+// 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/emcal_barrel_pions.cxx b/benchmarks/ecal/scripts/emcal_barrel_pions.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..02537b5603ed5912a956cbd961df3f93bdda0b7b
--- /dev/null
+++ b/benchmarks/ecal/scripts/emcal_barrel_pions.cxx
@@ -0,0 +1,80 @@
+//////////////////////////////////////////////////////////////
+// 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/emcal_barrel_pions_analysis.cxx b/benchmarks/ecal/scripts/emcal_barrel_pions_analysis.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c5c6f4421173a8bf169d8fcc8cb011835300072d
--- /dev/null
+++ b/benchmarks/ecal/scripts/emcal_barrel_pions_analysis.cxx
@@ -0,0 +1,171 @@
+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/CalorimeterHitCollection.h"
+#include "dd4pod/TrackerHitCollection.h"
+#include "eicd/RawCalorimeterHitCollection.h"
+#include "eicd/RawCalorimeterHitData.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 "TF1.h"
+#include "TH1D.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, {"RecoEcalBarrelAstroPixHitsXY"})
+ .Define("ErecZ", ErecZ, {"RecoEcalBarrelAstroPixHitsZ"})
+ .Define("ncluster", ncluster, {"EcalBarrelAstroPixClusters"})
+ .Define("Ecluster", Ecluster, {"EcalBarrelAstroPixClusters"})
+ .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/emal_pions_ErecZ.png");
+ c3->SaveAs("results/emal_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);
+ 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/emcal_barrel_pions_reader.cxx b/benchmarks/ecal/scripts/emcal_barrel_pions_reader.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..e03771482f146313dedab1c5d4391d80e751392a
--- /dev/null
+++ b/benchmarks/ecal/scripts/emcal_barrel_pions_reader.cxx
@@ -0,0 +1,125 @@
+//////////////////////////
+// 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/tracking/central_electrons.sh b/benchmarks/tracking/central_electrons.sh
index b355228554a1b19350d8da7b7f349ffc6383f796..77f77bd49fbf53d88cf8420fd82328f513f99ce2 100644
--- a/benchmarks/tracking/central_electrons.sh
+++ b/benchmarks/tracking/central_electrons.sh
@@ -95,7 +95,7 @@ fi
if [[ -z "${ANALYSIS_ONLY}" ]] ;
then
# Need to figure out how to pass file name to juggler from the commandline
- xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv gaudirun.py options/tracker_reconstruction.py
+ xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv gaudirun.py benchmarks/tracking/options/tracker_reconstruction.py
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running juggler"
exit 1
diff --git a/benchmarks/tracking/options/tracker_reconstruction.py b/benchmarks/tracking/options/tracker_reconstruction.py
index ede59d4f5d5e2f6b65af85f50d97ad734574623c..c587f4301ef751c99f7ef9f2eb19f613c311230c 100644
--- a/benchmarks/tracking/options/tracker_reconstruction.py
+++ b/benchmarks/tracking/options/tracker_reconstruction.py
@@ -46,7 +46,7 @@ from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
podioinput = PodioInput("PodioReader",
- collections=["mcparticles","SiTrackerEndcapHits","SiTrackerBarrelHits","EcalBarrelHits"])#, OutputLevel=DEBUG)"SiVertexBarrelHits",
+ collections=["mcparticles","SiTrackerEndcapHits","SiTrackerBarrelHits","EcalBarrelAstroPixHits"])#, OutputLevel=DEBUG)"SiVertexBarrelHits",
## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
copier = MCCopier("MCCopier",
@@ -57,8 +57,8 @@ trkcopier = TrkCopier("TrkCopier",
outputCollection="SiTrackerBarrelHits2")
ecal_digi = EMCalorimeterDigi("ecal_digi",
- inputHitCollection="EcalBarrelHits",
- outputHitCollection="RawEcalBarrelHits")
+ inputHitCollection="EcalBarrelAstroPixHits",
+ outputHitCollection="RawEcalBarrelAstroPixHits")
ufsd_digi = UFSDTrackerDigi("ufsd_digi",
inputHitCollection="SiTrackerBarrelHits",
@@ -76,13 +76,13 @@ ufsd_digi2 = UFSDTrackerDigi("ufsd_digi2",
ecal_reco = EMCalReconstruction("ecal_reco",
- inputHitCollection="RawEcalBarrelHits",
- outputHitCollection="RecEcalBarrelHits",
+ inputHitCollection="RawEcalBarrelAstroPixHits",
+ outputHitCollection="RecEcalBarrelAstroPixHits",
minModuleEdep=0.0*units.MeV,
OutputLevel=WARNING)
simple_cluster = SimpleClustering("simple_cluster",
- inputHitCollection="RecEcalBarrelHits",
+ inputHitCollection="RecEcalBarrelAstroPixHits",
outputClusters="SimpleClusters",
minModuleEdep=1.0*units.MeV,
maxDistance=50.0*units.cm,
diff --git a/options/tracker_reconstruction.py b/options/tracker_reconstruction.py
index 87130a7c491ea444ab5243e76c556b8e9d62c2cf..8fe85264697772791ed5e238527da7956f9da344 100644
--- a/options/tracker_reconstruction.py
+++ b/options/tracker_reconstruction.py
@@ -51,7 +51,7 @@ from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
podioinput = PodioInput("PodioReader",
- collections=["mcparticles","SiTrackerEndcapHits","SiTrackerBarrelHits","EcalBarrelHits"])#, OutputLevel=DEBUG)
+ collections=["mcparticles","SiTrackerEndcapHits","SiTrackerBarrelHits","EcalBarrelAstroPixHits"])#, OutputLevel=DEBUG)
#"SiVertexBarrelHits",
dummy = MC2DummyParticle("MC2Dummy",
@@ -67,8 +67,8 @@ trkcopier = TrkCopier("TrkCopier",
outputCollection="SiTrackerBarrelHits2")
ecal_digi = EMCalorimeterDigi("ecal_digi",
- inputHitCollection="EcalBarrelHits",
- outputHitCollection="RawEcalBarrelHits")
+ inputHitCollection="EcalBarrelAstroPixHits",
+ outputHitCollection="RawEcalBarrelAstroPixHits")
ufsd_digi = UFSDTrackerDigi("ufsd_digi",
inputHitCollection="SiTrackerBarrelHits",
@@ -88,13 +88,13 @@ ufsd_digi2 = UFSDTrackerDigi("ufsd_digi2",
ecal_reco = EMCalReconstruction("ecal_reco",
- inputHitCollection="RawEcalBarrelHits",
- outputHitCollection="RecEcalBarrelHits",
+ inputHitCollection="RawEcalBarrelAstroPixHits",
+ outputHitCollection="RecEcalBarrelAstroPixHits",
minModuleEdep=0.0*units.MeV)
#OutputLevel=DEBUG)
simple_cluster = SimpleClustering("simple_cluster",
- inputHitCollection="RecEcalBarrelHits",
+ inputHitCollection="RecEcalBarrelAstroPixHits",
outputClusters="SimpleClusters",
minModuleEdep=1.0*units.MeV,
maxDistance=50.0*units.cm)