diff --git a/benchmarks/clustering/options/full_cal_reco.py b/benchmarks/clustering/options/full_cal_reco.py
index 1c7d0508e61963abbd5f535a95c8e942c2b390e3..830726e4e113e8c2fbe2d90401706700788e0649 100644
--- a/benchmarks/clustering/options/full_cal_reco.py
+++ b/benchmarks/clustering/options/full_cal_reco.py
@@ -1,5 +1,5 @@
'''
- An example script to digitize/reconstruct calorimeter hits
+ An example option file to digitize/reconstruct/clustering calorimeter hits
'''
from Gaudi.Configuration import *
import os
@@ -14,7 +14,7 @@ detector_path = str(os.environ.get("JUGGLER_DETECTOR_PATH", "."))
compact_path = os.path.join(detector_path, detector_name)
# get sampling fractions from system environment variable, 1.0 by default
-ce_ecal_sf = float(os.environ.get("CE_ECAL_SAMP_FRAC", 0.253))
+ci_ecal_sf = float(os.environ.get("CI_ECAL_SAMP_FRAC", 0.253))
cb_ecal_sf = float(os.environ.get("CB_ECAL_SAMP_FRAC", 0.01324))
cb_hcal_sf = float(os.environ.get("CB_HCAL_SAMP_FRAC", 0.038))
ci_hcal_sf = float(os.environ.get("CI_HCAL_SAMP_FRAC", 0.025))
@@ -63,60 +63,62 @@ copier = MCCopier("MCCopier", inputCollection="mcparticles", outputCollection="m
# Crystal Endcap Ecal
+ce_ecal_daq = dict(
+ dynamicRangeADC=5.*GeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=3)
+
ce_ecal_digi = CalHitDigi("ce_ecal_digi",
inputHitCollection="EcalEndcapNHits",
outputHitCollection="EcalEndcapNHitsDigi",
energyResolutions=[0., 0.02, 0.],
- dynamicRangeADC=5.*GeV, # digi settings must match with reco
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=3)
+ **ce_ecal_daq)
ce_ecal_reco = CalHitReco("ce_ecal_reco",
inputHitCollection="EcalEndcapNHitsDigi",
outputHitCollection="EcalEndcapNHitsReco",
- dynamicRangeADC=5.*GeV,
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=3,
- thresholdFactor=4, # 4 sigma
- minimumHitEdep=1.0*MeV, # discard low energy hits
+ thresholdFactor=4, # 4 sigma cut on pedestal sigma
readoutClass="EcalEndcapNHits",
- sectorField="sector")
+ sectorField="sector",
+ **ce_ecal_daq)
ce_ecal_cl = IslandCluster("ce_ecal_cl",
# OutputLevel=DEBUG,
inputHitCollection="EcalEndcapNHitsReco",
outputHitCollection="EcalEndcapNClusterHits",
splitCluster=False,
+ minClusterHitEdep=1.0*MeV, # discard low energy hits
minClusterCenterEdep=30*MeV,
- groupRanges=[2.2*cm, 2.2*cm])
+ dimScaledDist=1.8) # dimension scaled dist is good for hybrid sectors with different module size
ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
- inputHitCollection="EcalEndcapNClusterHits",
- outputClusterCollection="EcalEndcapNClusters",
+ inputHitCollection="EcalEndcapNClusterHits",
+ outputClusterCollection="EcalEndcapNClusters",
samplingFraction=0.998, # this accounts for a small fraction of leakage
logWeightBase=4.6)
# Endcap Sampling Ecal
+ci_ecal_daq = dict(
+ dynamicRangeADC=50.*MeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=10)
+
ci_ecal_digi = CalHitDigi("ci_ecal_digi",
- inputHitCollection="EcalEndcapPHits",
- outputHitCollection="EcalEndcapPHitsDigi",
- dynamicRangeADC=50.*MeV,
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=10)
+ inputHitCollection="EcalEndcapPHits",
+ outputHitCollection="EcalEndcapPHitsDigi",
+ **ci_ecal_daq)
ci_ecal_reco = CalHitReco("ci_ecal_reco",
inputHitCollection="EcalEndcapPHitsDigi",
outputHitCollection="EcalEndcapPHitsReco",
- dynamicRangeADC=50.*MeV,
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=10,
- thresholdFactor=5.0)
+ thresholdFactor=5.0,
+ **ci_ecal_daq)
+
# merge hits in different layer (projection to local x-y plane)
ci_ecal_merger = CalHitsMerger("ci_ecal_merger",
+ # OutputLevel=DEBUG,
inputHitCollection="EcalEndcapPHitsReco",
outputHitCollection="EcalEndcapPHitsRecoXY",
fields=["layer", "slice"],
@@ -124,45 +126,49 @@ ci_ecal_merger = CalHitsMerger("ci_ecal_merger",
readoutClass="EcalEndcapPHits")
ci_ecal_cl = IslandCluster("ci_ecal_cl",
+ # OutputLevel=DEBUG,
inputHitCollection="EcalEndcapPHitsRecoXY",
outputHitCollection="EcalEndcapPClusterHits",
splitCluster=False,
- minClusterCenterEdep=30.*MeV,
- groupRanges=[5*mm, 5*mm])
+ minClusterCenterEdep=10.*MeV,
+ localDistXY=[10*mm, 10*mm])
ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
inputHitCollection="EcalEndcapPClusterHits",
- outputClusterCollection="EcalEndcapPClusters",
+ outputClusterCollection="EcalEndcapPClusters",
logWeightBase=6.2,
- samplingFraction=ce_ecal_sf)
+ samplingFraction=ci_ecal_sf)
# Central Barrel Ecal (Imaging Cal.)
-cb_ecal_digi = CalHitDigi("cb_ecal_digi",
- inputHitCollection="EcalBarrelHits",
- outputHitCollection="EcalBarrelHitsDigi",
- energyResolutions=[0., 0.02, 0.], # 2% flat resolution
+cb_ecal_daq = dict(
dynamicRangeADC=3*MeV,
capacityADC=8192,
pedestalMean=400,
pedestalSigma=20) # about 6 keV
+
+cb_ecal_digi = CalHitDigi("cb_ecal_digi",
+ inputHitCollection="EcalBarrelHits",
+ outputHitCollection="EcalBarrelHitsDigi",
+ energyResolutions=[0., 0.02, 0.], # 2% flat resolution
+ **cb_ecal_daq)
+
cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
inputHitCollection="EcalBarrelHitsDigi",
outputHitCollection="EcalBarrelHitsReco",
- dynamicRangeADC=3.*MeV,
- capacityADC=8192,
- pedestalMean=400,
- pedestalSigma=20,
thresholdFactor=3, # about 20 keV
readoutClass="EcalBarrelHits", # readout class
layerField="layer", # field to get layer id
- sectorField="module") # field to get sector id
+ sectorField="module", # field to get sector id
+ **cb_ecal_daq)
+
cb_ecal_cl = ImagingCluster("cb_ecal_cl",
inputHitCollection="EcalBarrelHitsReco",
outputHitCollection="EcalBarrelClusterHits",
- localRanges=[2.*mm, 2*mm], # same layer
- adjLayerRanges=[10*mrad, 10*mrad], # adjacent layer
- adjLayerDiff=2, # id diff for adjacent layer
- adjSectorDist=3.*cm) # different sector
+ localDistXY=[2.*mm, 2*mm], # same layer
+ layerDistEtaPhi=[10*mrad, 10*mrad], # adjacent layer
+ neighbourLayersRange=2, # id diff for adjacent layer
+ sectorDist=3.*cm) # different sector
+
cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
samplingFraction=cb_ecal_sf,
inputHitCollection="EcalBarrelClusterHits",
@@ -170,25 +176,26 @@ cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
outputLayerCollection="EcalBarrelLayers")
# Central Barrel Hcal
-cb_hcal_digi = CalHitDigi("cb_hcal_digi",
- inputHitCollection="HcalBarrelHits",
- outputHitCollection="HcalBarrelHitsDigi",
+cb_hcal_daq = dict(
dynamicRangeADC=50.*MeV,
capacityADC=32768,
pedestalMean=400,
pedestalSigma=10)
+cb_hcal_digi = CalHitDigi("cb_hcal_digi",
+ inputHitCollection="HcalBarrelHits",
+ outputHitCollection="HcalBarrelHitsDigi",
+ **cb_hcal_daq)
+
cb_hcal_reco = CalHitReco("cb_hcal_reco",
inputHitCollection="HcalBarrelHitsDigi",
outputHitCollection="HcalBarrelHitsReco",
- dynamicRangeADC=50.*MeV,
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=10,
thresholdFactor=5.0,
readoutClass="HcalBarrelHits",
layerField="layer",
- sectorField="module")
+ sectorField="module",
+ **cb_hcal_daq)
+
cb_hcal_merger = CalHitsMerger("cb_hcal_merger",
inputHitCollection="HcalBarrelHitsReco",
outputHitCollection="HcalBarrelHitsRecoXY",
@@ -201,32 +208,32 @@ cb_hcal_cl = IslandCluster("cb_hcal_cl",
outputHitCollection="HcalBarrelClusterHits",
splitCluster=False,
minClusterCenterEdep=30.*MeV,
- groupRanges=[15.*cm, 15.*cm])
+ localDistXY=[15.*cm, 15.*cm])
cb_hcal_clreco = RecoCoG("cb_hcal_clreco",
inputHitCollection="HcalBarrelClusterHits",
- outputClusterCollection="HcalBarrelClusters",
+ outputClusterCollection="HcalBarrelClusters",
logWeightBase=6.2,
samplingFraction=cb_hcal_sf)
# Hcal Hadron Endcap
-ci_hcal_digi = CalHitDigi("ci_hcal_digi",
- inputHitCollection="HcalHadronEndcapHits",
- outputHitCollection="HcalHadronEndcapHitsDigi",
+ci_hcal_daq = dict(
dynamicRangeADC=50.*MeV,
capacityADC=32768,
pedestalMean=400,
pedestalSigma=10)
+ci_hcal_digi = CalHitDigi("ci_hcal_digi",
+ inputHitCollection="HcalHadronEndcapHits",
+ outputHitCollection="HcalHadronEndcapHitsDigi",
+ **ci_hcal_daq)
ci_hcal_reco = CalHitReco("ci_hcal_reco",
inputHitCollection="HcalHadronEndcapHitsDigi",
outputHitCollection="HcalHadronEndcapHitsReco",
- dynamicRangeADC=50.*MeV,
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=10,
- thresholdFactor=5.0)
+ thresholdFactor=5.0,
+ **ci_hcal_daq)
+
ci_hcal_merger = CalHitsMerger("ci_hcal_merger",
inputHitCollection="HcalHadronEndcapHitsReco",
outputHitCollection="HcalHadronEndcapHitsRecoXY",
@@ -239,31 +246,32 @@ ci_hcal_cl = IslandCluster("ci_hcal_cl",
outputHitCollection="HcalHadronEndcapClusterHits",
splitCluster=False,
minClusterCenterEdep=30.*MeV,
- groupRanges=[15.*cm, 15.*cm])
+ localDistXY=[15.*cm, 15.*cm])
ci_hcal_clreco = RecoCoG("ci_hcal_clreco",
inputHitCollection="HcalHadronEndcapClusterHits",
- outputClusterCollection="HcalHadronEndcapClusters",
+ outputClusterCollection="HcalHadronEndcapClusters",
logWeightBase=6.2,
samplingFraction=ci_hcal_sf)
# Hcal Electron Endcap
-ce_hcal_digi = CalHitDigi("ce_hcal_digi",
- inputHitCollection="HcalElectronEndcapHits",
- outputHitCollection="HcalElectronEndcapHitsDigi",
+ce_hcal_daq = dict(
dynamicRangeADC=50.*MeV,
capacityADC=32768,
pedestalMean=400,
pedestalSigma=10)
+ce_hcal_digi = CalHitDigi("ce_hcal_digi",
+ inputHitCollection="HcalElectronEndcapHits",
+ outputHitCollection="HcalElectronEndcapHitsDigi",
+ **ce_hcal_daq)
+
ce_hcal_reco = CalHitReco("ce_hcal_reco",
inputHitCollection="HcalElectronEndcapHitsDigi",
outputHitCollection="HcalElectronEndcapHitsReco",
- dynamicRangeADC=50.*MeV,
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=10,
- thresholdFactor=5.0)
+ thresholdFactor=5.0,
+ **ce_hcal_daq)
+
ce_hcal_merger = CalHitsMerger("ce_hcal_merger",
inputHitCollection="HcalElectronEndcapHitsReco",
outputHitCollection="HcalElectronEndcapHitsRecoXY",
@@ -276,34 +284,34 @@ ce_hcal_cl = IslandCluster("ce_hcal_cl",
outputHitCollection="HcalElectronEndcapClusterHits",
splitCluster=False,
minClusterCenterEdep=30.*MeV,
- groupRanges=[15.*cm, 15.*cm])
+ localDistXY=[15.*cm, 15.*cm])
ce_hcal_clreco = RecoCoG("ce_hcal_clreco",
inputHitCollection="HcalElectronEndcapClusterHits",
- outputClusterCollection="HcalElectronEndcapClusters",
+ outputClusterCollection="HcalElectronEndcapClusters",
logWeightBase=6.2,
samplingFraction=ce_hcal_sf)
-podout.outputCommands = ['drop *', 'keep mcparticles*', 'keep *Reco', 'keep *Digi', 'keep *Clusters', 'keep *Layers', 'drop *ProtoClusters']
+podout.outputCommands = ['drop *',
+ 'keep mcparticles2',
+ 'keep *Digi',
+ 'keep *Reco*',
+ 'keep *ClusterHits',
+ 'keep *Clusters',
+ 'keep *Layers']
ApplicationMgr(
TopAlg = [podin, copier,
- ce_ecal_digi, ce_ecal_reco,
- ci_ecal_digi, ci_ecal_reco,
- cb_ecal_digi, cb_ecal_reco,
- cb_hcal_digi, cb_hcal_reco,
- ce_hcal_digi, ce_hcal_reco,
- ci_hcal_digi, ci_hcal_reco,
- ce_ecal_cl, ce_ecal_clreco,
- ci_ecal_merger, ci_ecal_cl, ci_ecal_clreco,
- cb_ecal_cl, cb_ecal_clreco,
- cb_hcal_merger, cb_hcal_cl, cb_hcal_clreco,
- ce_hcal_merger, ce_hcal_cl, ce_hcal_clreco,
- ci_hcal_merger, ci_hcal_cl, ci_hcal_clreco,
+ ce_ecal_digi, ce_ecal_reco, ce_ecal_cl, ce_ecal_clreco,
+ ci_ecal_digi, ci_ecal_reco, ci_ecal_merger, ci_ecal_cl, ci_ecal_clreco,
+ cb_ecal_digi, cb_ecal_reco, cb_ecal_cl, cb_ecal_clreco,
+ cb_hcal_digi, cb_hcal_reco, cb_hcal_merger, cb_hcal_cl, cb_hcal_clreco,
+ ce_hcal_digi, ce_hcal_reco, ce_hcal_merger, ce_hcal_cl, ce_hcal_clreco,
+ ci_hcal_digi, ci_hcal_reco, ci_hcal_merger, ci_hcal_cl, ci_hcal_clreco,
podout],
EvtSel = 'NONE',
- EvtMax = n_events,
+ EvtMax = n_events,
ExtSvc = [podioevent],
OutputLevel=DEBUG
)
diff --git a/benchmarks/clustering/scripts/cluster_plots.py b/benchmarks/clustering/scripts/cluster_plots.py
index 9d0bd63353adc50144b1f032fd7677fdc8138f12..6f70fa7ff88cfc803324e50e258a8418fa21198a 100644
--- a/benchmarks/clustering/scripts/cluster_plots.py
+++ b/benchmarks/clustering/scripts/cluster_plots.py
@@ -97,7 +97,7 @@ def general_clusters_figure(df, collection, save, min_nhits=3):
fig, axs = plt.subplots(2, 2, figsize=(16, 12), dpi=120)
labels = [
('Number of Hits', 'Counts'),
- (r'$E$ (MeV)', 'Counts'),
+ (r'$E$ (GeV)', 'Counts'),
(r'$\theta$ (rad)', 'Counts'),
(r'$\phi$ (rad)', 'Counts'),
]
diff --git a/benchmarks/clustering/scripts/gen_central_electrons.cxx b/benchmarks/clustering/scripts/deprecated/gen_central_electrons.cxx
similarity index 100%
rename from benchmarks/clustering/scripts/gen_central_electrons.cxx
rename to benchmarks/clustering/scripts/deprecated/gen_central_electrons.cxx
diff --git a/benchmarks/clustering/scripts/gen_central_pions.cxx b/benchmarks/clustering/scripts/deprecated/gen_central_pions.cxx
similarity index 100%
rename from benchmarks/clustering/scripts/gen_central_pions.cxx
rename to benchmarks/clustering/scripts/deprecated/gen_central_pions.cxx
diff --git a/benchmarks/ecal/config.yml b/benchmarks/ecal/config.yml
index 0632b5cba1e104dbf951f7b53b2e627aebdf4be2..c14cb9f66867ad0a3ef88a826d483006747484e4 100644
--- a/benchmarks/ecal/config.yml
+++ b/benchmarks/ecal/config.yml
@@ -1,42 +1,56 @@
-emcal_crystal_electrons:
+ecal_endcap_e_electrons:
extends: .rec_benchmark
stage: process
timeout: 8 hours
script:
- - bash benchmarks/ecal/emcal_electrons.sh
+ - python benchmarks/ecal/run_emcal_benchmarks.py endcap_e --particle "electron" -n 100 -t "endcap_e_electron"
+ --pmin 5.0 --pmax 5.0 --amin 135 --amax 175
-emcal_crystal_pi0s:
+ecal_endcap_i_electrons:
extends: .rec_benchmark
- stage: run
- timeout: 12 hours
+ stage: process
+ timeout: 8 hours
script:
- - bash benchmarks/ecal/emcal_pi0s.sh
+ - python benchmarks/ecal/run_emcal_benchmarks.py endcap_i --particle "electron" -n 100 -t "endcap_i_electron"
+ --pmin 5.0 --pmax 5.0 --amin 5 --amax 45
-emcal_barrel_electrons:
+ecal_barrel_electrons:
extends: .rec_benchmark
- timeout: 48 hours
- stage: run
+ stage: process
+ timeout: 8 hours
script:
- - bash benchmarks/ecal/run_emcal_barrel_electrons.sh
+ - python benchmarks/ecal/run_emcal_benchmarks.py barrel --particle "electron" -n 100 -t "barrel_electron"
+ --pmin 5.0 --pmax 5.0 --amin 45 --amax 135
-emcal_barrel_pions:
+ecal_endcap_e_pions:
extends: .rec_benchmark
- timeout: 48 hours
- stage: run
+ stage: process
+ timeout: 8 hours
script:
- - bash benchmarks/ecal/run_emcal_barrel_pions.sh
+ - python benchmarks/ecal/run_emcal_benchmarks.py endcap_e --particle "pion-" -n 100 -t "endcap_e_pion"
+ --pmin 5.0 --pmax 5.0 --amin 135 --amax 175
-full_emcal_barrel_electrons:
+ecal_endcap_i_pions:
extends: .rec_benchmark
- timeout: 48 hours
- stage: run
+ stage: process
+ timeout: 8 hours
+ script:
+ - python benchmarks/ecal/run_emcal_benchmarks.py endcap_i --particle "pion-" -n 100 -t "endcap_i_pion"
+ --pmin 5.0 --pmax 5.0 --amin 5 --amax 45
+
+ecal_barrel_pions:
+ extends: .rec_benchmark
+ stage: process
+ timeout: 8 hours
script:
- - bash benchmarks/ecal/full_emcal_electrons.sh
+ - python benchmarks/ecal/run_emcal_benchmarks.py barrel --particle "pion-" -n 100 -t "barrel_pion"
+ --pmin 5.0 --pmax 5.0 --amin 45 --amax 135
ecal_collect:
stage: collect
needs:
- - ["emcal_crystal_electrons", "emcal_crystal_pi0s", "emcal_barrel_electrons", "emcal_barrel_pions", "full_emcal_barrel_electrons"]
+ - ["ecal_endcap_e_electrons", "ecal_endcap_e_pions", "ecal_endcap_i_electrons", "ecal_endcap_i_pions",
+ "ecal_barrel_electrons", "ecal_barrel_pions"]
script:
- echo "Done collecting artifacts."
diff --git a/benchmarks/ecal/emcal_electrons.sh b/benchmarks/ecal/emcal_electrons.sh
deleted file mode 100644
index 67e082318e3ba4356b2f74b26f89afa5229dcd27..0000000000000000000000000000000000000000
--- a/benchmarks/ecal/emcal_electrons.sh
+++ /dev/null
@@ -1,116 +0,0 @@
-#!/bin/bash
-
-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
-if [[ ! -n "${JUGGLER_DETECTOR_PATH}" ]] ; then
- export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH}
-fi
-
-if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
- export JUGGLER_N_EVENTS=100
-fi
-
-if [[ ! -n "${E_start}" ]] ; then
- export E_start=0.0
-fi
-
-if [[ ! -n "${E_end}" ]] ; then
- export E_end=30.0
-fi
-
-export JUGGLER_FILE_NAME_TAG="emcal_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
-# note datasets is now only used to develop datasets.
-#git clone https://eicweb.phy.anl.gov/EIC/datasets.git datasets
-root -b -q "benchmarks/ecal/scripts/emcal_electrons.cxx(${JUGGLER_N_EVENTS}, ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running script"
- exit 1
-fi
-root -b -q "benchmarks/ecal/scripts/emcal_electrons_reader.cxx(${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running script"
- exit 1
-fi
-
-# run geant4 simulations
-npsim --runType batch \
- --part.minimalKineticEnergy 1000*GeV \
- -v WARNING \
- --numberOfEvents ${JUGGLER_N_EVENTS} \
- --compactFile ${JUGGLER_DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
- --inputFiles ${JUGGLER_FILE_NAME_TAG}.hepmc \
- --outputFile ${JUGGLER_SIM_FILE}
-# Need to figure out how to pass file name to juggler from the commandline
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running npdet"
- exit 1
-fi
-xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
- gaudirun.py benchmarks/ecal/options/crystal_calorimeter_reco.py
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running juggler"
- exit 1
-fi
-
-
-mkdir -p results
-
-#rootls ${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}
-root -b -q "benchmarks/ecal/scripts/rec_emcal_electrons_reader.C(${E_start}, ${E_end}, \"${JUGGLER_REC_FILE}\")"
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running analysis script"
- exit 1
-fi
-#root -b -q "ecal/scripts/makeplot.C(${E_start}, ${E_end}, \"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\", \"results/rec_${JUGGLER_FILE_NAME_TAG}.txt\")"
-#root -b -q "ecal/scripts/makeplot_input.C(\"${JUGGLER_DETECTOR}/${JUGGLER_SIM_FILE}\", \"results/sim_${JUGGLER_FILE_NAME_TAG}.txt\")"
-
-root -b -q "benchmarks/ecal/scripts/crystal_cal_electrons.cxx(\"${JUGGLER_REC_FILE}\")"
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running analysis script"
- exit 1
-fi
-
-root -b -q "benchmarks/ecal/scripts/emcal_electrons_analysis.cxx(\"${JUGGLER_REC_FILE}\")"
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running analysis script"
- exit 1
-fi
-
-# Script to generate energy resolution plots
-root -b -q "benchmarks/ecal/scripts/rec_emcal_resolution_analysis.cxx(\"${JUGGLER_REC_FILE}\")"
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running analysis script"
- exit 1
-fi
-
-#paste results/sim_${JUGGLER_FILE_NAME_TAG}.txt results/rec_${JUGGLER_FILE_NAME_TAG}.txt > results/eng_${JUGGLER_FILE_NAME_TAG}.txt
-#root -b -q "ecal/scripts/read_eng.C(\"results/eng_${JUGGLER_FILE_NAME_TAG}.root\", \"results/eng_${JUGGLER_FILE_NAME_TAG}.txt\")"
-#root -b -q "ecal/scripts/cal_eng_res.C(\"results/eng_${JUGGLER_FILE_NAME_TAG}.root\")"
-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/emcal_pi0s.sh b/benchmarks/ecal/emcal_pi0s.sh
deleted file mode 100644
index c3895aedd6ea6680d3a395ac815d4aafbb4d258e..0000000000000000000000000000000000000000
--- a/benchmarks/ecal/emcal_pi0s.sh
+++ /dev/null
@@ -1,86 +0,0 @@
-#!/bin/bash
-# Based on emcal_electrons.sh script
-
-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 read .local/bin/env.sh for more in-depth explanations of the variables
-## and how they can be controlled.
-
-if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
- export JUGGLER_N_EVENTS=100
-fi
-
-if [[ ! -n "${JUGGLER_DETECTOR_PATH}" ]] ; then
- export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH}
-fi
-
-# File names
-export JUGGLER_FILE_NAME_TAG="emcal_uniform_pi0s"
-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}"
-# Datasets
-#git clone https://eicweb.phy.anl.gov/EIC/datasets.git datasets
-root -b -q "benchmarks/ecal/scripts/emcal_pi0.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running analysis script"
- exit 1
-fi
-root -b -q "benchmarks/ecal/scripts/emcal_pi0_reader.cxx(\"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running analysis script"
- exit 1
-fi
-
-echo "CHECK POINT FOR GEANT4 SIMULATION"
-
-echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
-echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
-echo "JUGGLER_FILE_NAME_TAG = ${JUGGLER_FILE_NAME_TAG}"
-echo "JUGGLER_SIM_FILE = ${JUGGLER_SIM_FILE}"
-# run geant4 simulations
-npsim --runType batch \
- --part.minimalKineticEnergy 1000*GeV \
- -v WARNING \
- --numberOfEvents ${JUGGLER_N_EVENTS} \
- --compactFile ${JUGGLER_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
-
-# Need to figure out how to pass file name to juggler from the commandline
-gaudirun.py benchmarks/ecal/options/crystal_calorimeter_reco.py
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running juggler"
- exit 1
-fi
-
-
-mkdir -p results
-root -b -q "benchmarks/ecal/scripts/emcal_pion_anlaysis.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/full_emcal_electrons.sh b/benchmarks/ecal/full_emcal_electrons.sh
deleted file mode 100644
index b634d19e3441206eeab780ab9abc479d5be6bbbb..0000000000000000000000000000000000000000
--- a/benchmarks/ecal/full_emcal_electrons.sh
+++ /dev/null
@@ -1,82 +0,0 @@
-#!/bin/bash
-
-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
-if [[ ! -n "${JUGGLER_DETECTOR_PATH}" ]] ; then
- export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH}
-fi
-
-if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
- export JUGGLER_N_EVENTS=100
-fi
-
-if [[ ! -n "${E_start}" ]] ; then
- export E_start=0.5
-fi
-
-if [[ ! -n "${E_end}" ]] ; then
- export E_end=5.0
-fi
-
-export JUGGLER_FILE_NAME_TAG="emcal_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
-# note datasets is now only used to develop datasets.
-#git clone https://eicweb.phy.anl.gov/EIC/datasets.git datasets
-root -b -q "benchmarks/ecal/scripts/full_emcal_electrons.cxx(${JUGGLER_N_EVENTS}, ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running script"
- exit 1
-fi
-
-# run geant4 simulations
-npsim --runType batch \
- -v WARNING \
- --part.minimalKineticEnergy 1*GeV \
- --numberOfEvents ${JUGGLER_N_EVENTS} \
- --compactFile ${JUGGLER_DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
- --inputFiles ${JUGGLER_FILE_NAME_TAG}.hepmc \
- --outputFile ${JUGGLER_SIM_FILE}
-# Need to figure out how to pass file name to juggler from the commandline
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running npdet"
- exit 1
-fi
-
-xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
- gaudirun.py benchmarks/ecal/options/full_em_calorimeter_reco.py
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running juggler"
- exit 1
-fi
-
-root -b -q "benchmarks/ecal/scripts/full_emcal_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/options/barrel.py b/benchmarks/ecal/options/barrel.py
new file mode 100644
index 0000000000000000000000000000000000000000..f44bf58fca4c4e76ce84125f76d0b10269b3e6de
--- /dev/null
+++ b/benchmarks/ecal/options/barrel.py
@@ -0,0 +1,99 @@
+'''
+ An example script to digitize/reconstruct/clustering endcap ecal hits
+'''
+from Gaudi.Configuration import *
+import os
+import ROOT
+
+from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
+from Configurables import ApplicationMgr, EICDataSvc, PodioInput, PodioOutput, GeoSvc
+from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad
+
+detector_name = str(os.environ.get("JUGGLER_DETECTOR", "athena"))
+detector_path = str(os.environ.get("JUGGLER_DETECTOR_PATH", "."))
+compact_path = str(os.environ.get("JUGGLER_COMPACT_PATH", "{}.xml".format(os.path.join(detector_path, detector_name))))
+
+# input and output
+input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()]
+output_rec = str(os.environ["JUGGLER_REC_FILE"])
+n_events = int(os.environ["JUGGLER_N_EVENTS"])
+cb_ecal_sf = float(os.environ.get("CB_ECAL_SAMP_FRAC", 0.01324))
+
+# geometry service
+geo_service = GeoSvc("GeoSvc", detectors=[compact_path], OutputLevel=INFO)
+# data service
+podioevent = EICDataSvc("EventDataSvc", inputs=input_sims)
+
+
+# juggler components
+from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
+from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
+from Configurables import Jug__Reco__ImagingPixelReco as ImCalPixelReco
+from Configurables import Jug__Reco__ImagingTopoCluster as ImagingCluster
+from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
+
+# branches needed from simulation root file
+sim_coll = [
+ "mcparticles",
+ "EcalBarrelHits",
+]
+
+# input and output
+podin = PodioInput("PodioReader", collections=sim_coll)
+podout = PodioOutput("out", filename=output_rec)
+# copier needed to get around input --> output copy bug. So truth (mcparticles) can be saved in output file
+copier = MCCopier("MCCopier", inputCollection="mcparticles", outputCollection="mcparticles2")
+
+
+# Central Barrel Ecal (Imaging Cal.)
+cb_ecal_daq = dict(
+ dynamicRangeADC=3*MeV,
+ capacityADC=8192,
+ pedestalMean=400,
+ pedestalSigma=20) # about 6 keV
+
+cb_ecal_digi = CalHitDigi("cb_ecal_digi",
+ inputHitCollection="EcalBarrelHits",
+ outputHitCollection="EcalBarrelHitsDigi",
+ energyResolutions=[0., 0.02, 0.], # 2% flat resolution
+ **cb_ecal_daq)
+
+cb_ecal_reco = ImCalPixelReco("cb_ecal_reco",
+ inputHitCollection="EcalBarrelHitsDigi",
+ outputHitCollection="EcalBarrelHitsReco",
+ thresholdFactor=3, # about 20 keV
+ readoutClass="EcalBarrelHits", # readout class
+ layerField="layer", # field to get layer id
+ sectorField="module", # field to get sector id
+ **cb_ecal_daq)
+
+cb_ecal_cl = ImagingCluster("cb_ecal_cl",
+ inputHitCollection="EcalBarrelHitsReco",
+ outputHitCollection="EcalBarrelClusterHits",
+ localDistXY=[2.*mm, 2*mm], # same layer
+ layerDistEtaPhi=[10*mrad, 10*mrad], # adjacent layer
+ neighbourLayersRange=2, # id diff for adjacent layer
+ sectorDist=3.*cm) # different sector
+
+cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
+ samplingFraction=cb_ecal_sf,
+ inputHitCollection="EcalBarrelClusterHits",
+ outputClusterCollection="EcalBarrelClusters",
+ outputLayerCollection="EcalBarrelLayers")
+
+podout.outputCommands = ['drop *',
+ 'keep mcparticles2',
+ 'keep *HitsReco',
+ 'keep *HitsDigi',
+ 'keep *ClusterHits',
+ 'keep *Clusters']
+
+ApplicationMgr(
+ TopAlg = [podin, copier,
+ cb_ecal_digi, cb_ecal_reco, cb_ecal_cl, cb_ecal_clreco,
+ podout],
+ EvtSel = 'NONE',
+ EvtMax = n_events,
+ ExtSvc = [podioevent],
+ OutputLevel=DEBUG
+)
diff --git a/benchmarks/ecal/options/crystal_calorimeter_reco.py b/benchmarks/ecal/options/deprecated/crystal_calorimeter_reco.py
similarity index 100%
rename from benchmarks/ecal/options/crystal_calorimeter_reco.py
rename to benchmarks/ecal/options/deprecated/crystal_calorimeter_reco.py
diff --git a/benchmarks/ecal/options/emcal_barrel_reco.py b/benchmarks/ecal/options/deprecated/emcal_barrel_reco.py
similarity index 100%
rename from benchmarks/ecal/options/emcal_barrel_reco.py
rename to benchmarks/ecal/options/deprecated/emcal_barrel_reco.py
diff --git a/benchmarks/ecal/options/full_em_calorimeter_reco.py b/benchmarks/ecal/options/deprecated/full_em_calorimeter_reco.py
similarity index 100%
rename from benchmarks/ecal/options/full_em_calorimeter_reco.py
rename to benchmarks/ecal/options/deprecated/full_em_calorimeter_reco.py
diff --git a/benchmarks/ecal/options/endcap_e.py b/benchmarks/ecal/options/endcap_e.py
new file mode 100644
index 0000000000000000000000000000000000000000..9dee7fb8c6fa8a94643d8a030b1c6020b8433384
--- /dev/null
+++ b/benchmarks/ecal/options/endcap_e.py
@@ -0,0 +1,98 @@
+'''
+ An example script to digitize/reconstruct/clustering endcap ecal hits
+'''
+from Gaudi.Configuration import *
+import os
+import ROOT
+
+from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
+from Configurables import ApplicationMgr, EICDataSvc, PodioInput, PodioOutput, GeoSvc
+from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad
+
+detector_name = str(os.environ.get("JUGGLER_DETECTOR", "athena"))
+detector_path = str(os.environ.get("JUGGLER_DETECTOR_PATH", "."))
+compact_path = str(os.environ.get("JUGGLER_COMPACT_PATH", "{}.xml".format(os.path.join(detector_path, detector_name))))
+
+# input and output
+input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()]
+output_rec = str(os.environ["JUGGLER_REC_FILE"])
+n_events = int(os.environ["JUGGLER_N_EVENTS"])
+
+# geometry service
+geo_service = GeoSvc("GeoSvc", detectors=[compact_path], OutputLevel=INFO)
+# data service
+podioevent = EICDataSvc("EventDataSvc", inputs=input_sims)
+
+
+# juggler components
+from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
+from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
+from Configurables import Jug__Reco__CalorimeterHitReco as CalHitReco
+from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
+from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
+
+# branches needed from simulation root file
+sim_coll = [
+ "mcparticles",
+ "EcalEndcapNHits",
+]
+
+# input and output
+podin = PodioInput("PodioReader", collections=sim_coll)
+podout = PodioOutput("out", filename=output_rec)
+# copier needed to get around input --> output copy bug. So truth (mcparticles) can be saved in output file
+copier = MCCopier("MCCopier", inputCollection="mcparticles", outputCollection="mcparticles2")
+
+
+# Crystal Endcap Ecal
+ce_ecal_daq = dict(
+ dynamicRangeADC=5.*GeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=3)
+
+ce_ecal_digi = CalHitDigi("ce_ecal_digi",
+ inputHitCollection="EcalEndcapNHits",
+ outputHitCollection="EcalEndcapNHitsDigi",
+ energyResolutions=[0., 0.02, 0.],
+ **ce_ecal_daq)
+
+ce_ecal_reco = CalHitReco("ce_ecal_reco",
+ inputHitCollection="EcalEndcapNHitsDigi",
+ outputHitCollection="EcalEndcapNHitsReco",
+ thresholdFactor=4, # 4 sigma cut on pedestal sigma
+ readoutClass="EcalEndcapNHits",
+ sectorField="sector",
+ **ce_ecal_daq)
+
+ce_ecal_cl = IslandCluster("ce_ecal_cl",
+ # OutputLevel=DEBUG,
+ inputHitCollection="EcalEndcapNHitsReco",
+ outputHitCollection="EcalEndcapNClusterHits",
+ splitCluster=False,
+ minClusterHitEdep=1.0*MeV, # discard low energy hits
+ minClusterCenterEdep=30*MeV,
+ dimScaledDist=1.8) # hybrid calorimeter with different dimensions, using a scaled dist
+
+ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
+ inputHitCollection="EcalEndcapNClusterHits",
+ outputClusterCollection="EcalEndcapNClusters",
+ samplingFraction=0.998, # this accounts for a small fraction of leakage
+ logWeightBase=4.6)
+
+podout.outputCommands = ['drop *',
+ 'keep mcparticles2',
+ 'keep *HitsReco',
+ 'keep *HitsDigi',
+ 'keep *ClusterHits',
+ 'keep *Clusters']
+
+ApplicationMgr(
+ TopAlg = [podin, copier,
+ ce_ecal_digi, ce_ecal_reco, ce_ecal_cl, ce_ecal_clreco,
+ podout],
+ EvtSel = 'NONE',
+ EvtMax = n_events,
+ ExtSvc = [podioevent],
+ OutputLevel=DEBUG
+)
diff --git a/benchmarks/ecal/options/endcap_i.py b/benchmarks/ecal/options/endcap_i.py
new file mode 100644
index 0000000000000000000000000000000000000000..7e5553b609cd27cc1c00ba20c31501fbff9fc87c
--- /dev/null
+++ b/benchmarks/ecal/options/endcap_i.py
@@ -0,0 +1,106 @@
+
+'''
+ An example script to digitize/reconstruct/clustering endcap ecal hits
+'''
+from Gaudi.Configuration import *
+import os
+import ROOT
+
+from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
+from Configurables import ApplicationMgr, EICDataSvc, PodioInput, PodioOutput, GeoSvc
+from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad
+
+detector_name = str(os.environ.get("JUGGLER_DETECTOR", "athena"))
+detector_path = str(os.environ.get("JUGGLER_DETECTOR_PATH", "."))
+compact_path = str(os.environ.get("JUGGLER_COMPACT_PATH", "{}.xml".format(os.path.join(detector_path, detector_name))))
+ci_ecal_sf = float(os.environ.get("CI_ECAL_SAMP_FRAC", 0.253))
+
+# input and output
+input_sims = [f.strip() for f in str.split(os.environ["JUGGLER_SIM_FILE"], ",") if f.strip()]
+output_rec = str(os.environ["JUGGLER_REC_FILE"])
+n_events = int(os.environ["JUGGLER_N_EVENTS"])
+
+# geometry service
+geo_service = GeoSvc("GeoSvc", detectors=[compact_path], OutputLevel=INFO)
+# data service
+podioevent = EICDataSvc("EventDataSvc", inputs=input_sims)
+
+
+# juggler components
+from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
+from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
+from Configurables import Jug__Reco__CalorimeterHitReco as CalHitReco
+from Configurables import Jug__Reco__CalorimeterHitsMerger as CalHitsMerger
+from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
+from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
+
+# branches needed from simulation root file
+sim_coll = [
+ "mcparticles",
+ "EcalEndcapPHits",
+]
+
+# input and output
+podin = PodioInput("PodioReader", collections=sim_coll)
+podout = PodioOutput("out", filename=output_rec)
+# copier needed to get around input --> output copy bug. So truth (mcparticles) can be saved in output file
+copier = MCCopier("MCCopier", inputCollection="mcparticles", outputCollection="mcparticles2")
+
+
+# Endcap sampling Ecal
+ci_ecal_daq = dict(
+ dynamicRangeADC=50.*MeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=10)
+
+ci_ecal_digi = CalHitDigi("ci_ecal_digi",
+ inputHitCollection="EcalEndcapPHits",
+ outputHitCollection="EcalEndcapPHitsDigi",
+ **ci_ecal_daq)
+
+ci_ecal_reco = CalHitReco("ci_ecal_reco",
+ inputHitCollection="EcalEndcapPHitsDigi",
+ outputHitCollection="EcalEndcapPHitsReco",
+ thresholdFactor=5.0,
+ **ci_ecal_daq)
+
+# merge hits in different layer (projection to local x-y plane)
+ci_ecal_merger = CalHitsMerger("ci_ecal_merger",
+ # OutputLevel=DEBUG,
+ inputHitCollection="EcalEndcapPHitsReco",
+ outputHitCollection="EcalEndcapPHitsRecoXY",
+ fields=["layer", "slice"],
+ fieldRefNumbers=[1, 0],
+ readoutClass="EcalEndcapPHits")
+
+ci_ecal_cl = IslandCluster("ci_ecal_cl",
+ # OutputLevel=DEBUG,
+ inputHitCollection="EcalEndcapPHitsRecoXY",
+ outputHitCollection="EcalEndcapPClusterHits",
+ splitCluster=False,
+ minClusterCenterEdep=10.*MeV,
+ localDistXY=[10*mm, 10*mm])
+
+ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
+ inputHitCollection="EcalEndcapPClusterHits",
+ outputClusterCollection="EcalEndcapPClusters",
+ logWeightBase=6.2,
+ samplingFraction=ci_ecal_sf)
+
+podout.outputCommands = ['drop *',
+ 'keep mcparticles2',
+ 'keep *HitsReco*',
+ 'keep *HitsDigi',
+ 'keep *ClusterHits',
+ 'keep *Clusters']
+
+ApplicationMgr(
+ TopAlg = [podin, copier,
+ ci_ecal_digi, ci_ecal_reco, ci_ecal_merger, ci_ecal_cl, ci_ecal_clreco,
+ podout],
+ EvtSel = 'NONE',
+ EvtMax = n_events,
+ ExtSvc = [podioevent],
+ OutputLevel=DEBUG
+)
diff --git a/benchmarks/ecal/run_emcal_barrel_electrons.sh b/benchmarks/ecal/run_emcal_barrel_electrons.sh
deleted file mode 100755
index bd6351d8c9200ddb33c31e2a2ea068b170e7bd72..0000000000000000000000000000000000000000
--- a/benchmarks/ecal/run_emcal_barrel_electrons.sh
+++ /dev/null
@@ -1,98 +0,0 @@
-#!/bin/bash
-
-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
-if [[ ! -n "${JUGGLER_DETECTOR_PATH}" ]] ; then
- export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH}
-fi
-
-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
-
-if [[ ! -n "${CB_EMCAL_SAMP_FRAC}" ]] ; then
- export CB_EMCAL_SAMP_FRAC=0.014
-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 ${JUGGLER_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
-
-
-rootls -t "${JUGGLER_SIM_FILE}"
-
-# 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
deleted file mode 100755
index b25086a56bfcccabdb480173582b25de6d5cfb52..0000000000000000000000000000000000000000
--- a/benchmarks/ecal/run_emcal_barrel_pions.sh
+++ /dev/null
@@ -1,95 +0,0 @@
-#!/bin/bash
-
-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
-if [[ ! -n "${JUGGLER_DETECTOR_PATH}" ]] ; then
- export JUGGLER_DETECTOR_PATH=${DETECTOR_PATH}
-fi
-
-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
-
-if [[ ! -n "${CB_EMCAL_SAMP_FRAC}" ]] ; then
- export CB_EMCAL_SAMP_FRAC=0.01
-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 ${JUGGLER_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/run_emcal_benchmarks.py b/benchmarks/ecal/run_emcal_benchmarks.py
new file mode 100755
index 0000000000000000000000000000000000000000..c895735e6fb0917594803a1f3f8fc34dd7b205bd
--- /dev/null
+++ b/benchmarks/ecal/run_emcal_benchmarks.py
@@ -0,0 +1,121 @@
+#! /usr/local/bin/python3
+'''
+ A python script to help run emcal benchmarks
+ simulation -> reconstruction -> analysis
+ Author: Chao Peng, Jihee Kim (ANL)
+ Date: 07/04/2021
+'''
+import os
+import sys
+import subprocess
+import argparse
+
+
+# type: option files, analysis files
+default_type = {
+ 'endcap_e': [['endcap_e.py'], []],
+ 'endcap_i': [['endcap_i.py'], []],
+ 'barrel': [['barrel.py'], []],
+# 'all': [['all_ecal.py'], []],
+}
+
+default_compact = os.path.join(os.environ.get('JUGGLER_DETECTOR_PATH', os.environ.get('DETECTOR_PATH', '')),
+ '{}.xml'.format(os.environ.get('JUGGLER_DETECTOR', '')))
+parser = argparse.ArgumentParser()
+parser.add_argument('run_type', help='Run type, support {}'.format(list(default_type.keys())))
+parser.add_argument('-n', '--numberOfEvents', dest='nev', type=int, default=100, help='Number of events to process.')
+parser.add_argument('-t', '--nametag', type=str, default='IMCAL_ML', help='Name tag for output files.')
+parser.add_argument('--seed', type=int, default=-1, help='Random seed to scripts.')
+parser.add_argument('--process', type=str, default='sim, rec, ana', help='Processes to be executed (sim, rec, ana).')
+parser.add_argument('--numberOfLayers', dest='nlayer', type=int, default=20, help='number of layers in ML data.')
+parser.add_argument('--numberOfHits', dest='nhit', type=int, default=20, help='number of hits in ML data.')
+parser.add_argument('--particles', type=str, default='electron', help='Partcile names, separated by \",\".')
+parser.add_argument('--amin', type=float, default=5, help='Minimum polar angle of particles (degree).')
+parser.add_argument('--amax', type=float, default=175, help='Maximum polar angle of particles (degree).')
+parser.add_argument('--pmin', type=float, default=0.5, help='Minimum momentum of particles (GeV).')
+parser.add_argument('--pmax', type=float, default=10, help='Maximum momentum of particles (GeV).')
+parser.add_argument('--compact', type=str, default=default_compact, help='Path to detector compact file.')
+parser.add_argument('--uploadSizeLimit', type=float, default=10, help='Upper limit of file size (Mb) to be uploaded.')
+
+args = parser.parse_args()
+kwargs = vars(args)
+
+if args.run_type not in default_type:
+ print('ERROR unsupported run type {}, must in {}'.format(args.run_type, list(default_type.keys())))
+ exit(1)
+opt_scripts, ana_scripts = default_type[args.run_type]
+
+gen_file = 'gen_{nametag}_{pmin}_{pmax}.hepmc'.format(**kwargs)
+sim_file = 'sim_{nametag}_{pmin}_{pmax}.root'.format(**kwargs)
+rec_file = 'rec_{nametag}_{pmin}_{pmax}.root'.format(**kwargs)
+
+procs = [p.strip() for p in args.process.split(',')]
+sdir = os.path.dirname(os.path.realpath(__file__))
+
+if 'sim' in procs:
+ # generate particles
+ gen_cmd = ['python', os.path.join(sdir, 'scripts', 'gen_particles.py'), gen_file,
+ '-n', '{}'.format(args.nev),
+ '-s', '{}'.format(args.seed),
+ '--angmin', '{}'.format(args.amin), '--angmax', '{}'.format(args.amax),
+ '--pmin', '{}'.format(args.pmin), '--pmax', '{}'.format(args.pmax),
+ '--particles', args.particles]
+ subprocess.run(gen_cmd)
+ # simulation
+ sim_cmd = ['npsim',
+ '--part.minimalKineticEnergy', '1*TeV',
+ '--numberOfEvents', '{}'.format(args.nev),
+ '--runType', 'batch',
+ '--inputFiles', gen_file,
+ '--outputFile', sim_file,
+ '--compact', args.compact,
+ '-v', 'WARNING']
+ if args.seed > 0:
+ sim_cmd += ['--random.seed', args.seed]
+ print(sim_cmd)
+ return_code = subprocess.run(sim_cmd).returncode
+ if return_code is not None and return_code < 0:
+ print("ERROR running simulation!")
+ exit(1)
+ subprocess.run(['rootls', '-t', sim_file])
+
+
+if 'rec' in procs:
+ # export to environment variables (used to pass arguments to the option file)
+ os.environ['JUGGLER_SIM_FILE'] = sim_file
+ os.environ['JUGGLER_REC_FILE'] = rec_file
+ os.environ['JUGGLER_COMPACT_PATH'] = args.compact
+ os.environ['JUGGLER_N_EVENTS'] = '{}'.format(args.nev)
+
+ juggler_xenv = os.path.join(os.environ.get('JUGGLER_INSTALL_PREFIX', '../local'), 'Juggler.xenv')
+
+ for opt in opt_scripts:
+ rec_cmd = ['xenv', '-x', juggler_xenv, 'gaudirun.py', os.path.join(sdir, 'options', opt)]
+ return_code = subprocess.run(rec_cmd).returncode
+ if return_code is not None and return_code < 0:
+ print('ERROR running juggler ({})!'.format(opt))
+ exit(1)
+ process = subprocess.run(['rootls', '-t', rec_file])
+
+
+if 'ana' in procs:
+ os.makedirs('results', exist_ok=True)
+ for ana in ana_scripts:
+ ana_cmd = ['python', os.path.join(sdir, 'scripts', ana), rec_file, '-o', 'results']
+ return_code = subprocess.run(ana_cmd).returncode
+ if return_code is not None and return_code < 0:
+ print('ERROR running analysis ({})!'.format(ana))
+ exit(1)
+
+
+# upload generated data files if it was small enough (< 10 Mb)
+for upload in [rec_file]:
+ process = subprocess.run(['stat', '--format=%s', upload], stdout=subprocess.PIPE)
+ size = int(process.stdout)/1024./1024.
+ if size > args.uploadSizeLimit:
+ print('Skip uploading file \"{}\" because its size ({:.2f} Mb) > {:.2f} Mb'\
+ .format(upload, size, args.uploadSizeLimit))
+ else:
+ os.system('cp {} results/.'.format(upload))
+ print('Upload file \"{}\", size = {:.2f} Mb'.format(upload, size))
+
diff --git a/benchmarks/ecal/scripts/cal_eng_res.C b/benchmarks/ecal/scripts/deprecated/cal_eng_res.C
similarity index 100%
rename from benchmarks/ecal/scripts/cal_eng_res.C
rename to benchmarks/ecal/scripts/deprecated/cal_eng_res.C
diff --git a/benchmarks/ecal/scripts/crystal_cal_electrons.cxx b/benchmarks/ecal/scripts/deprecated/crystal_cal_electrons.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/crystal_cal_electrons.cxx
rename to benchmarks/ecal/scripts/deprecated/crystal_cal_electrons.cxx
diff --git a/benchmarks/ecal/scripts/emcal_barrel_electrons.cxx b/benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/emcal_barrel_electrons.cxx
rename to benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons.cxx
diff --git a/benchmarks/ecal/scripts/emcal_barrel_electrons_analysis.cxx b/benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons_analysis.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/emcal_barrel_electrons_analysis.cxx
rename to benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons_analysis.cxx
diff --git a/benchmarks/ecal/scripts/emcal_barrel_electrons_reader.cxx b/benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons_reader.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/emcal_barrel_electrons_reader.cxx
rename to benchmarks/ecal/scripts/deprecated/emcal_barrel_electrons_reader.cxx
diff --git a/benchmarks/ecal/scripts/emcal_barrel_pions.cxx b/benchmarks/ecal/scripts/deprecated/emcal_barrel_pions.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/emcal_barrel_pions.cxx
rename to benchmarks/ecal/scripts/deprecated/emcal_barrel_pions.cxx
diff --git a/benchmarks/ecal/scripts/emcal_barrel_pions_analysis.cxx b/benchmarks/ecal/scripts/deprecated/emcal_barrel_pions_analysis.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/emcal_barrel_pions_analysis.cxx
rename to benchmarks/ecal/scripts/deprecated/emcal_barrel_pions_analysis.cxx
diff --git a/benchmarks/ecal/scripts/emcal_barrel_pions_reader.cxx b/benchmarks/ecal/scripts/deprecated/emcal_barrel_pions_reader.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/emcal_barrel_pions_reader.cxx
rename to benchmarks/ecal/scripts/deprecated/emcal_barrel_pions_reader.cxx
diff --git a/benchmarks/ecal/scripts/emcal_electrons.cxx b/benchmarks/ecal/scripts/deprecated/emcal_electrons.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/emcal_electrons.cxx
rename to benchmarks/ecal/scripts/deprecated/emcal_electrons.cxx
diff --git a/benchmarks/ecal/scripts/emcal_electrons_analysis.cxx b/benchmarks/ecal/scripts/deprecated/emcal_electrons_analysis.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/emcal_electrons_analysis.cxx
rename to benchmarks/ecal/scripts/deprecated/emcal_electrons_analysis.cxx
diff --git a/benchmarks/ecal/scripts/emcal_electrons_reader.cxx b/benchmarks/ecal/scripts/deprecated/emcal_electrons_reader.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/emcal_electrons_reader.cxx
rename to benchmarks/ecal/scripts/deprecated/emcal_electrons_reader.cxx
diff --git a/benchmarks/ecal/scripts/emcal_pi0.cxx b/benchmarks/ecal/scripts/deprecated/emcal_pi0.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/emcal_pi0.cxx
rename to benchmarks/ecal/scripts/deprecated/emcal_pi0.cxx
diff --git a/benchmarks/ecal/scripts/emcal_pi0_reader.cxx b/benchmarks/ecal/scripts/deprecated/emcal_pi0_reader.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/emcal_pi0_reader.cxx
rename to benchmarks/ecal/scripts/deprecated/emcal_pi0_reader.cxx
diff --git a/benchmarks/ecal/scripts/emcal_pion_anlaysis.cxx b/benchmarks/ecal/scripts/deprecated/emcal_pion_anlaysis.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/emcal_pion_anlaysis.cxx
rename to benchmarks/ecal/scripts/deprecated/emcal_pion_anlaysis.cxx
diff --git a/benchmarks/ecal/scripts/full_emcal_electrons.cxx b/benchmarks/ecal/scripts/deprecated/full_emcal_electrons.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/full_emcal_electrons.cxx
rename to benchmarks/ecal/scripts/deprecated/full_emcal_electrons.cxx
diff --git a/benchmarks/ecal/scripts/full_emcal_electrons_analysis.cxx b/benchmarks/ecal/scripts/deprecated/full_emcal_electrons_analysis.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/full_emcal_electrons_analysis.cxx
rename to benchmarks/ecal/scripts/deprecated/full_emcal_electrons_analysis.cxx
diff --git a/benchmarks/ecal/scripts/makeplot.C b/benchmarks/ecal/scripts/deprecated/makeplot.C
similarity index 100%
rename from benchmarks/ecal/scripts/makeplot.C
rename to benchmarks/ecal/scripts/deprecated/makeplot.C
diff --git a/benchmarks/ecal/scripts/makeplot_input.C b/benchmarks/ecal/scripts/deprecated/makeplot_input.C
similarity index 100%
rename from benchmarks/ecal/scripts/makeplot_input.C
rename to benchmarks/ecal/scripts/deprecated/makeplot_input.C
diff --git a/benchmarks/ecal/scripts/makeplot_pi0.C b/benchmarks/ecal/scripts/deprecated/makeplot_pi0.C
similarity index 100%
rename from benchmarks/ecal/scripts/makeplot_pi0.C
rename to benchmarks/ecal/scripts/deprecated/makeplot_pi0.C
diff --git a/benchmarks/ecal/scripts/rdf_test.cxx b/benchmarks/ecal/scripts/deprecated/rdf_test.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/rdf_test.cxx
rename to benchmarks/ecal/scripts/deprecated/rdf_test.cxx
diff --git a/benchmarks/ecal/scripts/read_eng.C b/benchmarks/ecal/scripts/deprecated/read_eng.C
similarity index 100%
rename from benchmarks/ecal/scripts/read_eng.C
rename to benchmarks/ecal/scripts/deprecated/read_eng.C
diff --git a/benchmarks/ecal/scripts/rec_emcal_electrons_reader.C b/benchmarks/ecal/scripts/deprecated/rec_emcal_electrons_reader.C
similarity index 100%
rename from benchmarks/ecal/scripts/rec_emcal_electrons_reader.C
rename to benchmarks/ecal/scripts/deprecated/rec_emcal_electrons_reader.C
diff --git a/benchmarks/ecal/scripts/rec_emcal_resolution_analysis.cxx b/benchmarks/ecal/scripts/deprecated/rec_emcal_resolution_analysis.cxx
similarity index 100%
rename from benchmarks/ecal/scripts/rec_emcal_resolution_analysis.cxx
rename to benchmarks/ecal/scripts/deprecated/rec_emcal_resolution_analysis.cxx
diff --git a/benchmarks/ecal/scripts/gen_particles.py b/benchmarks/ecal/scripts/gen_particles.py
new file mode 100644
index 0000000000000000000000000000000000000000..2d0b81c8278d345cff479161f107f08f62c55bf1
--- /dev/null
+++ b/benchmarks/ecal/scripts/gen_particles.py
@@ -0,0 +1,104 @@
+import os
+from pyHepMC3 import HepMC3 as hm
+import numpy as np
+import argparse
+
+
+PARTICLES = {
+ "pion0": (111, 0.1349766), # pi0
+ "pion+": (211, 0.13957018), # pi+
+ "pion-": (-211, 0.13957018), # pi-
+ "kaon0": (311, 0.497648), # K0
+ "kaon+": (321, 0.493677), # K+
+ "kaon-": (-321, 0.493677), # K-
+ "proton": (2212, 0.938272), # proton
+ "neutron": (2112, 0.939565), # neutron
+ "electron": (11, 0.51099895e-3), # electron
+ "positron": (-11, 0.51099895e-3),# positron
+ "photon": (22, 0), # photon
+}
+
+
+def gen_event(p, theta, phi, pid, mass):
+ evt = hm.GenEvent(momentum_unit=hm.Units.MomentumUnit.GEV, length_unit=hm.Units.LengthUnit.MM)
+ # final state
+ state = 1
+ e0 = np.sqrt(p*p + mass*mass)
+ px = np.cos(phi)*np.sin(theta)
+ py = np.sin(phi)*np.sin(theta)
+ pz = np.cos(theta)
+
+ # beam
+ pbeam = hm.GenParticle(hm.FourVector(0, 0, 0, 0.938272), 2212, 4)
+ ebeam = hm.GenParticle(hm.FourVector(0, 0, e0, np.sqrt(e0*e0 + 0.511e-3*0.511e-3)), 11, 4)
+
+ # out particle
+ hout = hm.GenParticle(hm.FourVector(px*p, py*p, pz*p, e0), pid, state)
+
+ # vertex
+ vert = hm.GenVertex()
+ vert.add_particle_in(ebeam)
+ vert.add_particle_in(pbeam)
+ vert.add_particle_out(hout)
+ evt.add_vertex(vert)
+ return evt
+
+
+if __name__ == "__main__":
+ parser = argparse.ArgumentParser()
+
+ parser.add_argument('output', help='path to the output file')
+ parser.add_argument('-n', type=int, default=1000, dest='nev', help='number of events to generate')
+ parser.add_argument('-s', type=int, default=-1, dest='seed', help='seed for random generator')
+ parser.add_argument('--parray', type=str, default="", dest='parray',
+ help='an array of momenta in GeV, separated by \",\"')
+ parser.add_argument('--pmin', type=float, default=8.0, dest='pmin', help='minimum momentum in GeV')
+ parser.add_argument('--pmax', type=float, default=100.0, dest='pmax', help='maximum momentum in GeV')
+ parser.add_argument('--angmin', type=float, default=0.0, dest='angmin', help='minimum angle in degree')
+ parser.add_argument('--angmax', type=float, default=20.0, dest='angmax', help='maximum angle in degree')
+ parser.add_argument('--phmin', type=float, default=0.0, dest='phmin', help='minimum angle in degree')
+ parser.add_argument('--phmax', type=float, default=360.0, dest='phmax', help='maximum angle in degree')
+ parser.add_argument('--particles', type=str, default='electron', dest='particles',
+ help='particle names, support {}'.format(list(PARTICLES.keys())))
+
+ args = parser.parse_args()
+
+ # random seed (< 0 will get it from enviroment variable 'SEED', or a system random number)
+ if args.seed < 0:
+ args.seed = os.environ.get('SEED', int.from_bytes(os.urandom(4), byteorder='big', signed=False))
+ print("Random seed is {}".format(args.seed))
+ np.random.seed(args.seed)
+
+ output = hm.WriterAscii(args.output);
+ if output.failed():
+ print("Cannot open file \"{}\"".format(args.output))
+ sys.exit(2)
+
+ # build particle info
+ parts = []
+ for pid in args.particles.split(','):
+ pid = pid.strip()
+ if pid not in PARTICLES.keys():
+ print('pid {:d} not found in dictionary, ignored.'.format(pid))
+ continue
+ parts.append(PARTICLES[pid])
+
+ # p values
+ pvals = np.random.uniform(args.pmin, args.pmax, args.nev) if not args.parray else \
+ np.random.choice([float(p.strip()) for p in args.parray.split(',')], args.nev)
+ thvals = np.random.uniform(args.angmin, args.angmax, args.nev)/180.*np.pi
+ phivals = np.random.uniform(args.phmin, args.phmax, args.nev)/180.*np.pi
+ partvals = [parts[i] for i in np.random.choice(len(parts), args.nev)]
+
+ count = 0
+ for p, theta, phi, (pid, mass) in zip(pvals, thvals, phivals, partvals):
+ if (count % 1000 == 0):
+ print("Generated {} events".format(count), end='\r')
+ evt = gen_event(p, theta, phi, pid, mass)
+ output.write_event(evt)
+ evt.clear()
+ count += 1
+
+ print("Generated {} events".format(args.nev))
+ output.close()
+
diff --git a/benchmarks/imaging_ecal/options/imaging_topocluster.py b/benchmarks/imaging_ecal/options/imaging_topocluster.py
index 722985117bcb120eaca5921baf673751d766cb6b..fafb38156f47a6ebdaac9e3fc5ef308c0037a223 100644
--- a/benchmarks/imaging_ecal/options/imaging_topocluster.py
+++ b/benchmarks/imaging_ecal/options/imaging_topocluster.py
@@ -38,51 +38,58 @@ out = PodioOutput("out", filename=kwargs['output'])
podioinput = PodioInput("PodioReader", collections=["mcparticles", "EcalBarrelHits"], OutputLevel=DEBUG)
copier = MCCopier("MCCopier",
- inputCollection="mcparticles",
- outputCollection="mcparticles2",
- OutputLevel=DEBUG)
+ OutputLevel=DEBUG,
+ inputCollection="mcparticles",
+ outputCollection="mcparticles2")
calcopier = CalCopier("CalCopier",
- inputCollection="EcalBarrelHits",
- outputCollection="EcalBarrelHits2",
- OutputLevel=DEBUG)
+ OutputLevel=DEBUG,
+ inputCollection="EcalBarrelHits",
+ outputCollection="EcalBarrelHits2")
+
+# use the same daq_setting for digi/reco pair
+daq_setting = dict(
+ dynamicRangeADC=3*units.MeV,
+ capacityADC=32767,
+ pedestalMean=400,
+ pedestalSigma=50) # 50/32767*3 MeV ~ 5 keV
imcaldigi = CalorimeterHitDigi("imcal_digi",
- inputHitCollection="EcalBarrelHits",
- outputHitCollection="DigiEcalBarrelHits",
- inputEnergyUnit=units.GeV,
- inputTimeUnit=units.ns,
- energyResolutions=[0., 0.02, 0.],
- dynamicRangeADC=3*units.MeV,
- pedestalSigma=40,
- OutputLevel=DEBUG)
+ OutputLevel=DEBUG,
+ inputHitCollection="EcalBarrelHits",
+ outputHitCollection="DigiEcalBarrelHits",
+ energyResolutions=[0., 0.02, 0.],
+ **daq_setting)
imcalreco = ImagingPixelReco("imcal_reco",
- inputHitCollection="DigiEcalBarrelHits",
- outputHitCollection="RecoEcalBarrelHits",
- dynamicRangeADC=3.*units.MeV,
- pedestalSigma=40,
- readoutClass="EcalBarrelHits",
- layerField="layer",
- sectorField="module")
+ OutputLevel=DEBUG,
+ inputHitCollection="DigiEcalBarrelHits",
+ outputHitCollection="RecoEcalBarrelHits",
+ readoutClass="EcalBarrelHits",
+ layerField="layer",
+ sectorField="module",
+ **daq_setting)
+
imcalcluster = ImagingTopoCluster("imcal_cluster",
- inputHitCollection="RecoEcalBarrelHits",
- outputHitCollection="EcalBarrelClusterHits",
- localRanges=[2.*units.mm, 2*units.mm],
- adjLayerRanges=[10*units.mrad, 10*units.mrad],
- adjLayerDiff=2,
- adjSectorDist=3.*units.cm)
+ OutputLevel=DEBUG,
+ inputHitCollection="RecoEcalBarrelHits",
+ outputHitCollection="EcalBarrelClusterHits",
+ localDistXY=[2.*units.mm, 2*units.mm],
+ layerDistEtaPhi=[10*units.mrad, 10*units.mrad],
+ neighbourLayersRange=2,
+ sectorDist=3.*units.cm)
clusterreco = ImagingClusterReco("imcal_clreco",
- inputHitCollection="EcalBarrelClusterHits",
- outputLayerCollection="EcalBarrelClustersLayers",
- outputClusterCollection="EcalBarrelClusters",
- samplingFraction=sf,
- OutputLevel=DEBUG)
+ OutputLevel=DEBUG,
+ inputHitCollection="EcalBarrelClusterHits",
+ outputLayerCollection="EcalBarrelClustersLayers",
+ outputClusterCollection="EcalBarrelClusters",
+ samplingFraction=sf)
out.outputCommands = ["keep *"]
ApplicationMgr(
TopAlg=[podioinput, copier, calcopier, imcaldigi, imcalreco, imcalcluster, clusterreco, out],
EvtSel='NONE',
- EvtMax=kwargs['nev'],
+ # EvtMax=kwargs['nev'],
+ EvtMax=3,
ExtSvc=[podioevent],
OutputLevel=DEBUG
)
diff --git a/benchmarks/imaging_ecal/run_emcal_barrel.sh b/benchmarks/imaging_ecal/run_emcal_barrel.sh
index 43cc78c3fcf3523fb8da5241ee62cdf2f0b9b5f2..8584bef24e013b83a3a83808faae29ad43bbd15d 100644
--- a/benchmarks/imaging_ecal/run_emcal_barrel.sh
+++ b/benchmarks/imaging_ecal/run_emcal_barrel.sh
@@ -141,8 +141,6 @@ for iev in "${ADDR[@]}"; do
fi
done
-# TODO add this part back (solve mismatched EICD issue and add clusterreco back in option file)
-":" << COMMENT
python ${CB_EMCAL_SCRIPT_DIR}/energy_profile.py \
${CB_EMCAL_REC_FILE} --type=EM --energy=${CB_EMCAL_ENERGY} -o results/${particle} \
--save=results/profile.csv --color=royalblue
@@ -150,7 +148,6 @@ if [[ "$?" -ne "0" ]] ; then
echo "ERROR running analysis script: energy_profile"
exit 1
fi
-COMMENT
root_filesize=$(stat --format=%s "${CB_EMCAL_REC_FILE}")
if [[ "${CB_EMCAL_NUMEV}" -lt "500" ]] ; then
diff --git a/benchmarks/imaging_ecal/scripts/energy_profile.py b/benchmarks/imaging_ecal/scripts/energy_profile.py
index 234e646121657f6e958433f56f64aaa1a796cf62..d5b9b8aa21cd59ec2672abb6273f7d00fe6a57ea 100644
--- a/benchmarks/imaging_ecal/scripts/energy_profile.py
+++ b/benchmarks/imaging_ecal/scripts/energy_profile.py
@@ -47,7 +47,7 @@ if __name__ == '__main__':
dfe.loc[:, 'total_edep'] = dfe.groupby(['event', 'cluster'])['edep'].transform('sum')
# dfe.loc[:, 'efrac'] = dfe['edep']/dfe['total_edep']
dfe.loc[:, 'efrac'] = dfe['edep']/(args.energy*1000.)
- dfe = dfe[dfe['cluster'] == 1]
+ dfe = dfe[dfe['cluster'] == 0]
os.makedirs(args.outdir, exist_ok=True)
@@ -102,6 +102,8 @@ if __name__ == '__main__':
for ax, layer in zip(ax.flat, layers.flatten()):
data = dfe[dfe['layer'] == layer]
+ if not len(data):
+ continue
ax.hist(data['efrac'].values*100., weights=[1/float(len(data))]*len(data), bins=layers_bins,
ec='black', color=args.color)
ax.tick_params(labelsize=24)
diff --git a/benchmarks/imaging_ecal/scripts/utils.py b/benchmarks/imaging_ecal/scripts/utils.py
index ff8e368829eef58aa89e9585de3c05535a8e4ae0..8c9bfadfac14a56dcbf2b31ccef212088339b256 100644
--- a/benchmarks/imaging_ecal/scripts/utils.py
+++ b/benchmarks/imaging_ecal/scripts/utils.py
@@ -98,7 +98,7 @@ def get_hits_data(path, evnums=None, branch='RecoEcalBarrelHits'):
events.GetEntry(iev)
for hit in getattr(events, branch):
- dbuf[idb] = (iev, hit.clusterID, hit.layerID, hit.position.x, hit.position.y, hit.position.z, hit.edep)
+ dbuf[idb] = (iev, hit.clusterID, hit.layerID, hit.position.x, hit.position.y, hit.position.z, hit.edep*1000.)
idb += 1
return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'layer', 'x', 'y', 'z', 'edep'])
@@ -123,7 +123,7 @@ def get_layers_data(path, evnums=None, branch="EcalBarrelClustersLayers"):
events.GetEntry(iev)
for layer in getattr(events, branch):
dbuf[idb] = (iev, layer.clusterID, layer.layerID,
- layer.position.x, layer.position.y, layer.position.z, layer.edep)
+ layer.position.x, layer.position.y, layer.position.z, layer.edep*1000.)
idb += 1
return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'layer', 'x', 'y', 'z', 'edep'])
@@ -147,7 +147,7 @@ def get_clusters_data(path, evnums=None, branch='EcalBarrelClustersReco'):
events.GetEntry(iev)
for k, cl in enumerate(getattr(events, branch)):
- dbuf[idb] = (iev, k, cl.nhits, cl.edep, cl.cl_theta, cl.cl_phi)
+ dbuf[idb] = (iev, k, cl.nhits, cl.edep*1000., cl.cl_theta, cl.cl_phi)
idb += 1
return pd.DataFrame(data=dbuf[:idb], columns=['event', 'cluster', 'nhits', 'edep', 'cl_theta', 'cl_phi'])