From 88bbb45318a5486a986c6d93f8b4fbd6be4dea75 Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Fri, 18 Jun 2021 01:23:49 +0000
Subject: [PATCH] split clustering script to two: digi and clustering
---
benchmarks/clustering/full_cal_clusters.sh | 63 ++--
.../clustering/options/full_cal_clustering.py | 307 ------------------
.../clustering/options/full_cal_clusters.py | 201 ++++++++++++
.../clustering/options/full_cal_digi.py | 188 +++++++++++
.../clustering/scripts/cluster_plots.py | 24 +-
5 files changed, 440 insertions(+), 343 deletions(-)
delete mode 100644 benchmarks/clustering/options/full_cal_clustering.py
create mode 100644 benchmarks/clustering/options/full_cal_clusters.py
create mode 100644 benchmarks/clustering/options/full_cal_digi.py
diff --git a/benchmarks/clustering/full_cal_clusters.sh b/benchmarks/clustering/full_cal_clusters.sh
index 54498afd..220a6360 100644
--- a/benchmarks/clustering/full_cal_clusters.sh
+++ b/benchmarks/clustering/full_cal_clusters.sh
@@ -35,7 +35,7 @@ do
shift # past value
;;
-n|--nevents)
- export FULL_CAL_NUMEV="$2"
+ export JUGGLER_N_EVENTS="$2"
shift # past argument
shift # past value
;;
@@ -59,31 +59,32 @@ do
done
set -- "${POSITIONAL[@]}" # restore positional parameters
-export FULL_CAL_COMPACT_PATH=${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml
+export JUGGLER_COMPACT_PATH=${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml
-if [[ ! -n "${FULL_CAL_NUMEV}" ]] ; then
- export FULL_CAL_NUMEV=1000
+if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
+ export JUGGLER_N_EVENTS=1000
fi
if [[ ! -n "${FULL_CAL_PMIN}" ]] ; then
- export FULL_CAL_ENERGY=1.0
+ export FULL_CAL_PMIN=1.0
fi
if [[ ! -n "${FULL_CAL_PMAX}" ]] ; then
- export FULL_CAL_ENERGY=10.0
+ export FULL_CAL_PMIN=10.0
fi
-export FULL_CAL_NAME_TAG="${nametag}"
-export FULL_CAL_GEN_FILE="${FULL_CAL_NAME_TAG}.hepmc"
+export JUGGLER_FILE_NAME_TAG="${nametag}"
+export JUGGLER_GEN_FILE="gen_${JUGGLER_FILE_NAME_TAG}.hepmc"
-export FULL_CAL_SIM_FILE="sim_${FULL_CAL_NAME_TAG}.root"
-export FULL_CAL_REC_FILE="rec_${FULL_CAL_NAME_TAG}.root"
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_DIGI_FILE="digi_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
-echo "FULL_CAL_NUMEV = ${FULL_CAL_NUMEV}"
-echo "FULL_CAL_COMPACT_PATH = ${FULL_CAL_COMPACT_PATH}"
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "JUGGLER_COMPACT_PATH = ${JUGGLER_COMPACT_PATH}"
# Generate the input events
-python benchmarks/imaging_ecal/scripts/gen_particles.py ${FULL_CAL_GEN_FILE} -n ${FULL_CAL_NUMEV}\
+python benchmarks/imaging_ecal/scripts/gen_particles.py ${JUGGLER_GEN_FILE} -n ${JUGGLER_N_EVENTS}\
--angmin 5 --angmax 175 --phmin 0 --phmax 360 --pmin ${FULL_CAL_PMIN} --pmax ${FULL_CAL_PMAX}\
--particles="${particle}"
@@ -92,23 +93,23 @@ if [[ "$?" -ne "0" ]] ; then
exit 1
fi
-ls -lh ${FULL_CAL_GEN_FILE}
+ls -lh ${JUGGLER_GEN_FILE}
# Run geant4 simulations
npsim --runType batch \
-v WARNING \
--part.minimalKineticEnergy "1*TeV" \
- --numberOfEvents ${FULL_CAL_NUMEV} \
- --compactFile ${FULL_CAL_COMPACT_PATH} \
- --inputFiles ${FULL_CAL_GEN_FILE} \
- --outputFile ${FULL_CAL_SIM_FILE}
+ --numberOfEvents ${JUGGLER_N_EVENTS} \
+ --compactFile ${JUGGLER_COMPACT_PATH} \
+ --inputFiles ${JUGGLER_GEN_FILE} \
+ --outputFile ${JUGGLER_SIM_FILE}
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running npdet"
exit 1
fi
-rootls -t "${FULL_CAL_SIM_FILE}"
+rootls -t "${JUGGLER_SIM_FILE}"
# Directory for results
mkdir -p results
@@ -117,21 +118,31 @@ FULL_CAL_OPTION_DIR=benchmarks/clustering/options
FULL_CAL_SCRIPT_DIR=benchmarks/clustering/scripts
# Run Juggler
+# Digitization
xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
- gaudirun.py ${FULL_CAL_OPTION_DIR}/full_cal_clustering.py
+ gaudirun.py ${FULL_CAL_OPTION_DIR}/full_cal_digi.py
if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running juggler"
+ echo "ERROR running digitization (juggler)"
exit 1
fi
-# run analysis scripts
-python ${FULL_CAL_SCRIPT_DIR}/cluster_plots.py ${FULL_CAL_REC_FILE} -o results
+# Clustering
+xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+ gaudirun.py ${FULL_CAL_OPTION_DIR}/full_cal_clusters.py
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running clustering (juggler)"
+ exit 1
+fi
+
+
+# Run analysis scripts
+python ${FULL_CAL_SCRIPT_DIR}/cluster_plots.py ${JUGGLER_SIM_FILE} ${JUGGLER_REC_FILE} -o results
-root_filesize=$(stat --format=%s "${FULL_CAL_REC_FILE}")
-if [[ "${FULL_CAL_NUMEV}" -lt "500" ]] ; then
+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 ${FULL_CAL_REC_FILE} results/.
+ cp ${JUGGLER_REC_FILE} results/.
fi
fi
diff --git a/benchmarks/clustering/options/full_cal_clustering.py b/benchmarks/clustering/options/full_cal_clustering.py
deleted file mode 100644
index 13647aa7..00000000
--- a/benchmarks/clustering/options/full_cal_clustering.py
+++ /dev/null
@@ -1,307 +0,0 @@
-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("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))
-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))
-hcal_hadronendcap_sf = float(os.environ.get("CE_HCAL_HADRONENDCAP_SAMP_FRAC", 0.025))
-hcal_electronendcap_sf = float(os.environ.get("CE_HCAL_ELECTRONENDCAP_SAMP_FRAC", 0.025))
-
-# input and output
-input_sims = [f.strip() for f in str.split(os.environ["FULL_CAL_SIM_FILE"], ",") if f.strip()]
-output_rec = str(os.environ["FULL_CAL_REC_FILE"])
-n_events = int(os.environ["FULL_CAL_NUMEV"])
-
-# geometry service
-geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(compact_path)])
-# 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__ImagingPixelReco as ImCalPixelReco
-from Configurables import Jug__Reco__ImagingTopoCluster as ImagingCluster
-
-from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
-from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
-
-# branches needed from simulation root file
-sim_coll = [
- "mcparticles",
- "CrystalEcalHits",
- "EcalEndcapHits",
- "EcalBarrelHits",
- "HcalBarrelHits",
- "HcalHadronEndcapHits",
- "HcalElectronEndcapHits",
-]
-
-# 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_digi = CalHitDigi("ce_ecal_digi",
- inputHitCollection="CrystalEcalHits",
- outputHitCollection="CrystalEcalHitsDigi",
- energyResolutions=[0., 0.02, 0.],
- dynamicRangeADC=5.*GeV, # digi settings must match with reco
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=3)
-
-ce_ecal_reco = CalHitReco("ce_ecal_reco",
- inputHitCollection="CrystalEcalHitsDigi",
- outputHitCollection="CrystalEcalHitsReco",
- dynamicRangeADC=5.*GeV,
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=3,
- thresholdFactor=4, # 4 sigma
- minimumHitEdep=1.0*MeV) # discard low energy hits
-
-ce_ecal_cl = IslandCluster("ce_ecal_cl",
- # OutputLevel=DEBUG,
- inputHitCollection="CrystalEcalHitsReco",
- outputClusterCollection="CrystalEcalClusters",
- splitHitCollection="CrystalEcalHitsSplit",
- splitCluster=False,
- minClusterCenterEdep=30*MeV,
- groupRanges=[2.2*cm, 2.2*cm])
-
-ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
- clusterCollection="CrystalEcalClusters",
- samplingFraction=0.998, # this accounts for a small fraction of leakage
- logWeightBase=4.6)
-
-
-# Endcap Sampling Ecal
-ce_ecal2_digi = CalHitDigi("ce_ecal2_digi",
- inputHitCollection="EcalEndcapHits",
- outputHitCollection="EcalEndcapHitsDigi",
- dynamicRangeADC=50.*MeV,
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=10)
-
-ce_ecal2_reco = CalHitReco("ce_ecal2_reco",
- inputHitCollection="EcalEndcapHitsDigi",
- outputHitCollection="EcalEndcapHitsReco",
- dynamicRangeADC=50.*MeV,
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=10,
- thresholdFactor=5.0)
-
-# merge hits in different layer (projection to local x-y plane)
-ce_ecal2_merger = CalHitsMerger("ce_ecal2_merger",
- inputHitCollection="EcalEndcapHitsReco",
- outputHitCollection="EcalEndcapHitsRecoXY",
- fields=["layer", "slice"],
- fieldRefNumbers=[1, 0],
- readoutClass="EcalEndcapHits")
-
-ce_ecal2_cl = IslandCluster("ce_ecal2_cl",
- inputHitCollection="EcalEndcapHitsRecoXY",
- outputClusterCollection="EcalEndcapClusters",
- splitHitCollection="EcalEndcapHitsSplit",
- splitCluster=False,
- minClusterCenterEdep=30.*MeV,
- groupRanges=[5*mm, 5*mm])
-
-ce_ecal2_clreco = RecoCoG("ce_ecal2_clreco",
- clusterCollection="EcalEndcapClusters",
- logWeightBase=6.2,
- samplingFraction=ce_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
- dynamicRangeADC=3*MeV,
- capacityADC=8192,
- pedestalMean=400,
- pedestalSigma=20) # about 6 keV
-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
-cb_ecal_cl = ImagingCluster("cb_ecal_cl",
- inputHitCollection="EcalBarrelHitsReco",
- outputClusterCollection="EcalBarrelClusters",
- 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
-cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
- samplingFraction=cb_ecal_sf,
- inputClusterCollection="EcalBarrelClusters",
- outputLayerCollection="EcalBarrelLayers")
-
-
-# Central Barrel Hcal
-cb_hcal_digi = CalHitDigi("cb_hcal_digi",
- inputHitCollection="HcalBarrelHits",
- outputHitCollection="HcalBarrelHitsDigi",
- dynamicRangeADC=50.*MeV,
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=10)
-
-cb_hcal_reco = CalHitReco("cb_hcal_reco",
- inputHitCollection="HcalBarrelHitsDigi",
- outputHitCollection="HcalBarrelHitsReco",
- dynamicRangeADC=50.*MeV,
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=10,
- thresholdFactor=5.0)
-
-# merge hits in different layer (projection to local x-y plane)
-cb_hcal_merger = CalHitsMerger("cb_hcal_merger",
- inputHitCollection="HcalBarrelHitsReco",
- outputHitCollection="HcalBarrelHitsRecoXY",
- readoutClass="HcalBarrelHits",
- fields=["layer", "slice"],
- fieldRefNumbers=[1, 0])
-
-cb_hcal_cl = IslandCluster("cb_hcal_cl",
- inputHitCollection="HcalBarrelHitsRecoXY",
- outputClusterCollection="HcalBarrelClusters",
- splitHitCollection="HcalBarrelHitsSplit",
- splitCluster=False,
- minClusterCenterEdep=30.*MeV,
- groupRanges=[15.*cm, 15.*cm])
-
-cb_hcal_clreco = RecoCoG("cb_hcal_clreco",
- clusterCollection="HcalBarrelClusters",
- logWeightBase=6.2,
- samplingFraction=cb_hcal_sf)
-
-
-##Hadron Endcap HCal
-hcal_hadronendcap_digi = CalHitDigi("hcal_hadronendcap_digi",
- inputHitCollection="HcalHadronEndcapHits",
- outputHitCollection="HcalHadronEndcapHitsDigi",
- dynamicRangeADC=50.*MeV,
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=10)
-
-hcal_hadronendcap_reco = CalHitReco("hcal_hadronendcap_reco",
- inputHitCollection="HcalHadronEndcapHitsDigi",
- outputHitCollection="HcalHadronEndcapHitsReco",
- dynamicRangeADC=50.*MeV,
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=10,
- thresholdFactor=5.0)
-
-hcal_hadronendcap_merger = CalHitsMerger("hcal_hadronendcap_merger",
- inputHitCollection="HcalHadronEndcapHitsReco",
- outputHitCollection="HcalHadronEndcapHitsRecoXY",
- readoutClass="HcalHadronEndcapHits",
- fields=["layer", "slice"],
- fieldRefNumbers=[1, 0])
-
-hcal_hadronendcap_cl = IslandCluster("hcal_hadronendcap_cl",
- inputHitCollection="HcalHadronEndcapHitsRecoXY",
- outputClusterCollection="HcalHadronEndcapClusters",
- splitHitCollection="HcalHadronEndcapHitsSplit",
- splitCluster=False,
- minClusterCenterEdep=30.*MeV,
- groupRanges=[15.*cm, 15.*cm])
-
-hcal_hadronendcap_clreco = RecoCoG("hcal_hadronendcap_clreco",
- clusterCollection="HcalHadronEndcapClusters",
- logWeightBase=6.2,
- samplingFraction=hcal_hadronendcap_sf)
-
-#Hcal Electron Endcap
-hcal_electronendcap_digi = CalHitDigi("hcal_electronendcap_digi",
- inputHitCollection="HcalElectronEndcapHits",
- outputHitCollection="HcalElectronEndcapHitsDigi",
- dynamicRangeADC=50.*MeV,
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=10)
-
-hcal_electronendcap_reco = CalHitReco("hcal_electronendcap_reco",#
- inputHitCollection="HcalElectronEndcapHitsDigi",
- outputHitCollection="HcalElectronEndcapHitsReco",#
- dynamicRangeADC=50.*MeV,#
- capacityADC=32768,
- pedestalMean=400,
- pedestalSigma=10,
- thresholdFactor=5.0)
-hcal_electronendcap_merger = CalHitsMerger("hcal_electronendcap_merger",
- inputHitCollection="HcalElectronEndcapHitsReco",
- outputHitCollection="HcalElectronEndcapHitsRecoXY",
- readoutClass="HcalElectronEndcapHits",
- fields=["layer", "slice"],
- fieldRefNumbers=[1, 0])
-
-hcal_electronendcap_cl = IslandCluster("hcal_electronendcap_cl",
- inputHitCollection="HcalElectronEndcapHitsRecoXY",
- outputClusterCollection="HcalElectronEndcapClusters",
- splitHitCollection="HcalElectronEndcapHitsSplit",
- splitCluster=False,
- minClusterCenterEdep=30.*MeV,
- groupRanges=[15.*cm, 15.*cm])
-
-hcal_electronendcap_clreco = RecoCoG("hcal_electronendcap_clreco",
- clusterCollection="HcalElectronEndcapClusters",
- logWeightBase=6.2,
- samplingFraction=hcal_electronendcap_sf)
-
-
-podout.outputCommands = ["keep *"]
-
-ApplicationMgr(
- TopAlg = [podin, copier,
- ce_ecal_digi, ce_ecal_reco, ce_ecal_cl, ce_ecal_clreco,
- ce_ecal2_digi, ce_ecal2_reco, ce_ecal2_merger, ce_ecal2_cl, ce_ecal2_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,
- hcal_hadronendcap_digi, hcal_hadronendcap_reco, hcal_hadronendcap_merger,
- hcal_hadronendcap_cl, hcal_hadronendcap_clreco,
- hcal_electronendcap_digi,
- hcal_electronendcap_reco,
- hcal_electronendcap_merger,
- hcal_electronendcap_cl,
- #hcal_electronendcap_clreco,
- podout],
- EvtSel = 'NONE',
- EvtMax = n_events,
- ExtSvc = [podioevent],
- OutputLevel=DEBUG
-)
-
diff --git a/benchmarks/clustering/options/full_cal_clusters.py b/benchmarks/clustering/options/full_cal_clusters.py
new file mode 100644
index 00000000..c091624c
--- /dev/null
+++ b/benchmarks/clustering/options/full_cal_clusters.py
@@ -0,0 +1,201 @@
+'''
+ An example script to cluster reconstructed calorimeter 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("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))
+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))
+ce_hcal_sf = float(os.environ.get("CE_HCAL_SAMP_FRAC", 0.025))
+
+# input and output
+input_digis = [f.strip() for f in str.split(os.environ["JUGGLER_DIGI_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=["{}.xml".format(compact_path)])
+# data service
+podioevent = EICDataSvc("EventDataSvc", inputs=input_digis)
+
+
+# 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__ImagingPixelReco as ImCalPixelReco
+from Configurables import Jug__Reco__ImagingTopoCluster as ImagingCluster
+
+from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
+from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
+
+# branches needed from simulation root file
+digi_coll = [
+# "mcparticles2",
+ "CrystalEcalHitsReco",
+ "EcalEndcapHitsReco",
+ "EcalBarrelHitsReco",
+ "HcalBarrelHitsReco",
+ "HcalHadronEndcapHitsReco",
+ "HcalElectronEndcapHitsReco",
+]
+
+# input and output
+podin = PodioInput("PodioReader", collections=digi_coll, OutputLevel=DEBUG)
+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="mcparticles2", outputCollection="mcparticles3")
+
+
+# Crystal Endcap Ecal
+ce_ecal_cl = IslandCluster("ce_ecal_cl",
+ # OutputLevel=DEBUG,
+ inputHitCollection="CrystalEcalHitsReco",
+ outputClusterCollection="CrystalEcalClusters",
+ splitHitCollection="CrystalEcalHitsSplit",
+ splitCluster=False,
+ minClusterCenterEdep=30*MeV,
+ groupRanges=[2.2*cm, 2.2*cm])
+
+ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
+ clusterCollection="CrystalEcalClusters",
+ samplingFraction=0.998, # this accounts for a small fraction of leakage
+ logWeightBase=4.6)
+
+
+# Endcap Sampling Ecal
+# merge hits in different layer (projection to local x-y plane)
+ce_ecal2_merger = CalHitsMerger("ce_ecal2_merger",
+ inputHitCollection="EcalEndcapHitsReco",
+ outputHitCollection="EcalEndcapHitsRecoXY",
+ fields=["layer", "slice"],
+ fieldRefNumbers=[1, 0],
+ readoutClass="EcalEndcapHits")
+
+ce_ecal2_cl = IslandCluster("ce_ecal2_cl",
+ inputHitCollection="EcalEndcapHitsRecoXY",
+ outputClusterCollection="EcalEndcapClusters",
+ splitHitCollection="EcalEndcapHitsSplit",
+ splitCluster=False,
+ minClusterCenterEdep=30.*MeV,
+ groupRanges=[5*mm, 5*mm])
+
+ce_ecal2_clreco = RecoCoG("ce_ecal2_clreco",
+ clusterCollection="EcalEndcapClusters",
+ logWeightBase=6.2,
+ samplingFraction=ce_ecal_sf)
+
+
+# Central Barrel Ecal (Imaging Cal.)
+cb_ecal_cl = ImagingCluster("cb_ecal_cl",
+ inputHitCollection="EcalBarrelHitsReco",
+ outputClusterCollection="EcalBarrelClusters",
+ 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
+cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco",
+ samplingFraction=cb_ecal_sf,
+ inputClusterCollection="EcalBarrelClusters",
+ outputLayerCollection="EcalBarrelLayers")
+
+
+# Central Barrel Hcal
+# merge hits in different layer (projection to local x-y plane)
+cb_hcal_merger = CalHitsMerger("cb_hcal_merger",
+ inputHitCollection="HcalBarrelHitsReco",
+ outputHitCollection="HcalBarrelHitsRecoXY",
+ readoutClass="HcalBarrelHits",
+ fields=["layer", "slice"],
+ fieldRefNumbers=[1, 0])
+
+cb_hcal_cl = IslandCluster("cb_hcal_cl",
+ inputHitCollection="HcalBarrelHitsRecoXY",
+ outputClusterCollection="HcalBarrelClusters",
+ splitHitCollection="HcalBarrelHitsSplit",
+ splitCluster=False,
+ minClusterCenterEdep=30.*MeV,
+ groupRanges=[15.*cm, 15.*cm])
+
+cb_hcal_clreco = RecoCoG("cb_hcal_clreco",
+ clusterCollection="HcalBarrelClusters",
+ logWeightBase=6.2,
+ samplingFraction=cb_hcal_sf)
+
+
+# Hcal Hadron Endcap HCal
+ci_hcal_merger = CalHitsMerger("ci_hcal_merger",
+ inputHitCollection="HcalHadronEndcapHitsReco",
+ outputHitCollection="HcalHadronEndcapHitsRecoXY",
+ readoutClass="HcalHadronEndcapHits",
+ fields=["layer", "slice"],
+ fieldRefNumbers=[1, 0])
+
+ci_hcal_cl = IslandCluster("ci_hcal_cl",
+ inputHitCollection="HcalHadronEndcapHitsRecoXY",
+ outputClusterCollection="HcalHadronEndcapClusters",
+ splitHitCollection="HcalHadronEndcapHitsSplit",
+ splitCluster=False,
+ minClusterCenterEdep=30.*MeV,
+ groupRanges=[15.*cm, 15.*cm])
+
+ci_hcal_clreco = RecoCoG("ci_hcal_clreco",
+ clusterCollection="HcalHadronEndcapClusters",
+ logWeightBase=6.2,
+ samplingFraction=ci_hcal_sf)
+
+# Hcal Electron Endcap
+ce_hcal_merger = CalHitsMerger("ce_hcal_merger",
+ inputHitCollection="HcalElectronEndcapHitsReco",
+ outputHitCollection="HcalElectronEndcapHitsRecoXY",
+ readoutClass="HcalElectronEndcapHits",
+ fields=["layer", "slice"],
+ fieldRefNumbers=[1, 0])
+
+ce_hcal_cl = IslandCluster("ce_hcal_cl",
+ inputHitCollection="HcalElectronEndcapHitsRecoXY",
+ outputClusterCollection="HcalElectronEndcapClusters",
+ splitHitCollection="HcalElectronEndcapHitsSplit",
+ splitCluster=False,
+ minClusterCenterEdep=30.*MeV,
+ groupRanges=[15.*cm, 15.*cm])
+
+ce_hcal_clreco = RecoCoG("ce_hcal_clreco",
+ clusterCollection="HcalElectronEndcapClusters",
+ logWeightBase=6.2,
+ samplingFraction=ce_hcal_sf)
+
+
+podout.outputCommands = ["drop *", "keep mcparticles", "keep *Clusters", "keep *Layers"]
+
+ApplicationMgr(
+ TopAlg = [podin, # copier,
+ ce_ecal_cl, ce_ecal_clreco,
+ ce_ecal2_merger, ce_ecal2_cl, ce_ecal2_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,
+ podout],
+ EvtSel = 'NONE',
+ EvtMax = n_events,
+ ExtSvc = [podioevent],
+ OutputLevel=DEBUG
+)
+
diff --git a/benchmarks/clustering/options/full_cal_digi.py b/benchmarks/clustering/options/full_cal_digi.py
new file mode 100644
index 00000000..dcc584a1
--- /dev/null
+++ b/benchmarks/clustering/options/full_cal_digi.py
@@ -0,0 +1,188 @@
+'''
+ An example script to digitize/reconstruct calorimeter 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("DETECTOR_PATH", "."))
+compact_path = 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_DIGI_FILE"])
+n_events = int(os.environ["JUGGLER_N_EVENTS"])
+
+# geometry service
+geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(compact_path)])
+# 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__ImagingPixelReco as ImCalPixelReco
+from Configurables import Jug__Reco__ImagingTopoCluster as ImagingCluster
+
+from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
+from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
+
+# branches needed from simulation root file
+sim_coll = [
+ "mcparticles",
+ "CrystalEcalHits",
+ "EcalEndcapHits",
+ "EcalBarrelHits",
+ "HcalBarrelHits",
+ "HcalHadronEndcapHits",
+ "HcalElectronEndcapHits",
+]
+
+# 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_digi = CalHitDigi("ce_ecal_digi",
+ inputHitCollection="CrystalEcalHits",
+ outputHitCollection="CrystalEcalHitsDigi",
+ energyResolutions=[0., 0.02, 0.],
+ dynamicRangeADC=5.*GeV, # digi settings must match with reco
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=3)
+
+ce_ecal_reco = CalHitReco("ce_ecal_reco",
+ inputHitCollection="CrystalEcalHitsDigi",
+ outputHitCollection="CrystalEcalHitsReco",
+ dynamicRangeADC=5.*GeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=3,
+ thresholdFactor=4, # 4 sigma
+ minimumHitEdep=1.0*MeV) # discard low energy hits
+
+# Endcap Sampling Ecal
+ce_ecal2_digi = CalHitDigi("ce_ecal2_digi",
+ inputHitCollection="EcalEndcapHits",
+ outputHitCollection="EcalEndcapHitsDigi",
+ dynamicRangeADC=50.*MeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=10)
+
+ce_ecal2_reco = CalHitReco("ce_ecal2_reco",
+ inputHitCollection="EcalEndcapHitsDigi",
+ outputHitCollection="EcalEndcapHitsReco",
+ dynamicRangeADC=50.*MeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=10,
+ thresholdFactor=5.0)
+
+# Central Barrel Ecal (Imaging Cal.)
+cb_ecal_digi = CalHitDigi("cb_ecal_digi",
+ inputHitCollection="EcalBarrelHits",
+ outputHitCollection="EcalBarrelHitsDigi",
+ energyResolutions=[0., 0.02, 0.], # 2% flat resolution
+ dynamicRangeADC=3*MeV,
+ capacityADC=8192,
+ pedestalMean=400,
+ pedestalSigma=20) # about 6 keV
+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
+
+# Central Barrel Hcal
+cb_hcal_digi = CalHitDigi("cb_hcal_digi",
+ inputHitCollection="HcalBarrelHits",
+ outputHitCollection="HcalBarrelHitsDigi",
+ dynamicRangeADC=50.*MeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=10)
+
+cb_hcal_reco = CalHitReco("cb_hcal_reco",
+ inputHitCollection="HcalBarrelHitsDigi",
+ outputHitCollection="HcalBarrelHitsReco",
+ dynamicRangeADC=50.*MeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=10,
+ thresholdFactor=5.0)
+
+# Hcal Hadron Endcap
+ci_hcal_digi = CalHitDigi("ci_hcal_digi",
+ inputHitCollection="HcalHadronEndcapHits",
+ outputHitCollection="HcalHadronEndcapHitsDigi",
+ dynamicRangeADC=50.*MeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=10)
+
+ci_hcal_reco = CalHitReco("ci_hcal_reco",
+ inputHitCollection="HcalHadronEndcapHitsDigi",
+ outputHitCollection="HcalHadronEndcapHitsReco",
+ dynamicRangeADC=50.*MeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=10,
+ thresholdFactor=5.0)
+
+# Hcal Electron Endcap
+ce_hcal_digi = CalHitDigi("ce_hcal_digi",
+ inputHitCollection="HcalElectronEndcapHits",
+ outputHitCollection="HcalElectronEndcapHitsDigi",
+ dynamicRangeADC=50.*MeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=10)
+
+ce_hcal_reco = CalHitReco("ce_hcal_reco",
+ inputHitCollection="HcalElectronEndcapHitsDigi",
+ outputHitCollection="HcalElectronEndcapHitsReco",
+ dynamicRangeADC=50.*MeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=10,
+ thresholdFactor=5.0)
+
+
+podout.outputCommands = ['drop *', 'keep mcparticles2', 'keep *Reco', 'keep *Digi']
+
+ApplicationMgr(
+ TopAlg = [podin, copier,
+ ce_ecal_digi, ce_ecal_reco,
+ ce_ecal2_digi, ce_ecal2_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,
+ podout],
+ EvtSel = 'NONE',
+ EvtMax = n_events,
+ ExtSvc = [podioevent],
+ OutputLevel=DEBUG
+)
+
diff --git a/benchmarks/clustering/scripts/cluster_plots.py b/benchmarks/clustering/scripts/cluster_plots.py
index 71c59671..b7bcb95b 100644
--- a/benchmarks/clustering/scripts/cluster_plots.py
+++ b/benchmarks/clustering/scripts/cluster_plots.py
@@ -21,10 +21,10 @@ def load_root_macros(arg_macros):
print('Loading root macros: \'{}\' does not exist, skip loading it.'.format(path))
-def thrown_particles_figure(df, save):
+def thrown_particles_figure(df, save, mcbranch="mcparticles"):
# define truth particle info
- df = df.Define('isThrown', 'mcparticles2.genStatus == 1')\
- .Define('thrownParticles', 'mcparticles2[isThrown]')\
+ df = df.Define('isThrown', '{}.genStatus == 1'.format(mcbranch))\
+ .Define('thrownParticles', '{}[isThrown]'.format(mcbranch))\
.Define('thrownP', 'fourvec(thrownParticles)')\
.Define('thrownPID', 'pid(thrownParticles)')\
.Define('thrownEta', 'eta(thrownParticles)')\
@@ -116,10 +116,13 @@ def general_clusters_figure(df, collection, save, min_nhits=3):
if __name__ == '__main__':
parser = argparse.ArgumentParser()
- parser.add_argument('file', help='Path to reconstruction file.')
+ parser.add_argument('sim_file', help='Path to simulation output file.')
+ parser.add_argument('rec_file', help='Path to reconstruction output file.')
parser.add_argument('-o', dest='outdir', default='.', help='Output directory.')
parser.add_argument('-m', '--macros', dest='macros', default='rootlogon.C',
help='Macro files to be loaded by root, separated by \",\".')
+ parser.add_argument('--mc-branch', dest='mc', default='mcparticles',
+ help='Branch name of MC particles truth info.')
args = parser.parse_args()
# multi-threading for RDataFrame
@@ -136,15 +139,16 @@ if __name__ == '__main__':
# load macros (add libraries/headers root cannot automatically found here)
load_root_macros(args.macros)
- rdf = ROOT.RDataFrame('events', args.file)
+ rdf_sim = ROOT.RDataFrame('events', args.sim_file)
+ rdf_rec = ROOT.RDataFrame('events', args.rec_file)
- thrown_particles_figure(rdf, save=os.path.join(args.outdir, 'thrown_particles.png'))
- general_clusters_figure(rdf, collection='CrystalEcalClusters', min_nhits=5,
+ thrown_particles_figure(rdf_sim, save=os.path.join(args.outdir, 'thrown_particles.png'), mcbranch=args.mc)
+ general_clusters_figure(rdf_rec, collection='CrystalEcalClusters', min_nhits=5,
save=os.path.join(args.outdir, 'crystal_ecal_clusters.png'))
- general_clusters_figure(rdf, collection='EcalEndcapClusters', min_nhits=10,
+ general_clusters_figure(rdf_rec, collection='EcalEndcapClusters', min_nhits=10,
save=os.path.join(args.outdir, 'ecal_endcap_clusters.png'))
- general_clusters_figure(rdf, collection='EcalBarrelClusters',
+ general_clusters_figure(rdf_rec, collection='EcalBarrelClusters',
save=os.path.join(args.outdir, 'ecal_barrel_clusters.png'))
- general_clusters_figure(rdf, collection='HcalBarrelClusters',
+ general_clusters_figure(rdf_rec, collection='HcalBarrelClusters',
save=os.path.join(args.outdir, 'hcal_barrel_clusters.png'))
--
GitLab