From e6639cae5aea94f484ec1c7d01f5e7042b7dae4e Mon Sep 17 00:00:00 2001
From: Wouter Deconinck <wouter.deconinck@umanitoba.ca>
Date: Tue, 28 Sep 2021 17:44:17 +0000
Subject: [PATCH] ecal and hcal hits reconstruction
---
benchmarks/dis/dis.sh | 21 ++--
benchmarks/dvcs/dvcs.sh | 9 +-
benchmarks/dvmp/dvmp.sh | 22 ++---
benchmarks/u_omega/u_omega.sh | 9 +-
options/reconstruction.ecal.py | 168 ++++++++++++++++++++++++++++++++
options/reconstruction.hcal.py | 170 +++++++++++++++++++++++++++++++++
6 files changed, 370 insertions(+), 29 deletions(-)
create mode 100644 options/reconstruction.ecal.py
create mode 100644 options/reconstruction.hcal.py
diff --git a/benchmarks/dis/dis.sh b/benchmarks/dis/dis.sh
index a01af459..9ceb0de8 100755
--- a/benchmarks/dis/dis.sh
+++ b/benchmarks/dis/dis.sh
@@ -84,19 +84,16 @@ echo "Running the digitization and reconstruction"
## - JUGGLER_DETECTOR: detector package (part of global environment)
export JUGGLER_SIM_FILE=${SIM_FILE}
export JUGGLER_REC_FILE=${REC_FILE}
-xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
- gaudirun.py options/reconstruction.py
-## on-error, first retry running juggler again as there is still a random
-## crash we need to address FIXME
+for rec in options/*.py ; do
+ unset tag
+ [[ $(basename ${rec} .py) =~ (.*)\.(.*) ]] && tag=".${BASH_REMATCH[2]}"
+ JUGGLER_REC_FILE=${JUGGLER_REC_FILE/.root/${tag:-}.root} \
+ xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+ gaudirun.py ${rec}
+done
if [ "$?" -ne "0" ] ; then
- echo "Juggler crashed, retrying..."
- xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
- gaudirun.py options/reconstruction.py \
- 2>&1 > ${REC_LOG}
- if [ "$?" -ne "0" ] ; then
- echo "ERROR running juggler, both attempts failed"
- exit 1
- fi
+ echo "ERROR running juggler"
+ exit 1
fi
## =============================================================================
diff --git a/benchmarks/dvcs/dvcs.sh b/benchmarks/dvcs/dvcs.sh
index 1dd20afb..7a47ac14 100644
--- a/benchmarks/dvcs/dvcs.sh
+++ b/benchmarks/dvcs/dvcs.sh
@@ -115,8 +115,13 @@ fi
### Step 3. Run the reconstruction (juggler)
if [[ -n "${DO_REC}" || -n "${DO_ALL}" ]] ; then
- xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
- gaudirun.py options/reconstruction.py
+ for rec in options/*.py ; do
+ unset tag
+ [[ $(basename ${rec} .py) =~ (.*)\.(.*) ]] && tag=".${BASH_REMATCH[2]}"
+ JUGGLER_REC_FILE=${JUGGLER_REC_FILE/.root/${tag:-}.root} \
+ xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+ gaudirun.py ${rec}
+ done
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running juggler"
exit 1
diff --git a/benchmarks/dvmp/dvmp.sh b/benchmarks/dvmp/dvmp.sh
index 58a289fb..e9d608e1 100755
--- a/benchmarks/dvmp/dvmp.sh
+++ b/benchmarks/dvmp/dvmp.sh
@@ -84,20 +84,16 @@ echo "Running the digitization and reconstruction"
## - JUGGLER_DETECTOR: detector package (part of global environment)
export JUGGLER_SIM_FILE=${SIM_FILE}
export JUGGLER_REC_FILE=${REC_FILE}
-xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
- gaudirun.py options/reconstruction.py
-## 2>&1 > ${REC_LOG}
-## on-error, first retry running juggler again as there is still a random
-## crash we need to address FIXME
+for rec in options/*.py ; do
+ unset tag
+ [[ $(basename ${rec} .py) =~ (.*)\.(.*) ]] && tag=".${BASH_REMATCH[2]}"
+ JUGGLER_REC_FILE=${JUGGLER_REC_FILE/.root/${tag:-}.root} \
+ xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+ gaudirun.py ${rec}
+done
if [ "$?" -ne "0" ] ; then
- echo "Juggler crashed, retrying..."
- xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
- gaudirun.py options/reconstruction.py \
- 2>&1 > ${REC_LOG}
- if [ "$?" -ne "0" ] ; then
- echo "ERROR running juggler, both attempts failed"
- exit 1
- fi
+ echo "ERROR running juggler, both attempts failed"
+ exit 1
fi
## =============================================================================
diff --git a/benchmarks/u_omega/u_omega.sh b/benchmarks/u_omega/u_omega.sh
index fa0f69bb..6d051589 100644
--- a/benchmarks/u_omega/u_omega.sh
+++ b/benchmarks/u_omega/u_omega.sh
@@ -114,8 +114,13 @@ fi
### Step 3. Run the reconstruction (juggler)
if [[ -n "${DO_REC}" || -n "${DO_ALL}" ]] ; then
- xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
- gaudirun.py options/reconstruction.py
+ for rec in options/*.py ; do
+ unset tag
+ [[ $(basename ${rec} .py) =~ (.*)\.(.*) ]] && tag=".${BASH_REMATCH[2]}"
+ JUGGLER_REC_FILE=${JUGGLER_REC_FILE/.root/${tag:-}.root} \
+ xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+ gaudirun.py ${rec}
+ done
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running juggler"
exit 1
diff --git a/options/reconstruction.ecal.py b/options/reconstruction.ecal.py
new file mode 100644
index 00000000..9407fae8
--- /dev/null
+++ b/options/reconstruction.ecal.py
@@ -0,0 +1,168 @@
+from Gaudi.Configuration import *
+
+from Configurables import ApplicationMgr, AuditorSvc, EICDataSvc, PodioOutput, GeoSvc
+from GaudiKernel import SystemOfUnits as units
+from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad
+
+import json
+
+detector_name = "athena"
+if "JUGGLER_DETECTOR" in os.environ :
+ detector_name = str(os.environ["JUGGLER_DETECTOR"])
+
+detector_path = ""
+if "DETECTOR_PATH" in os.environ :
+ detector_path = str(os.environ["DETECTOR_PATH"])
+
+compact_path = os.path.join(detector_path, detector_name)
+
+# CAL reconstruction
+# get sampling fractions from system environment variable
+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"])
+
+# services
+services = []
+# geometry service
+services.append(GeoSvc("GeoSvc", detectors=["{}.xml".format(compact_path)], OutputLevel=WARNING))
+# data service
+services.append(EICDataSvc("EventDataSvc", inputs=input_sims, OutputLevel=WARNING))
+
+# juggler components
+from Configurables import PodioInput
+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",
+ "EcalEndcapNHits",
+ "EcalEndcapPHits",
+]
+
+# list of algorithms
+algorithms = []
+
+# input
+podin = PodioInput("PodioReader", collections=sim_coll)
+algorithms.append(podin)
+
+# Crystal Endcap Ecal
+ce_ecal_daq = dict(
+ dynamicRangeADC=5.*units.GeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=3)
+
+ce_ecal_digi = CalHitDigi("ce_ecal_digi",
+ inputHitCollection="EcalEndcapNHits",
+ outputHitCollection="EcalEndcapNRawHits",
+ energyResolutions=[0., 0.02, 0.],
+ **ce_ecal_daq)
+algorithms.append(ce_ecal_digi)
+
+ce_ecal_reco = CalHitReco("ce_ecal_reco",
+ inputHitCollection=ce_ecal_digi.outputHitCollection,
+ outputHitCollection="EcalEndcapNRecHits",
+ thresholdFactor=4, # 4 sigma cut on pedestal sigma
+ readoutClass="EcalEndcapNHits",
+ sectorField="sector",
+ **ce_ecal_daq)
+algorithms.append(ce_ecal_reco)
+
+ce_ecal_cl = IslandCluster("ce_ecal_cl",
+ inputHitCollection=ce_ecal_reco.outputHitCollection,
+ outputProtoClusterCollection="EcalEndcapNProtoClusters",
+ splitCluster=False,
+ minClusterHitEdep=1.0*units.MeV, # discard low energy hits
+ minClusterCenterEdep=30*units.MeV,
+ sectorDist=5.0*units.cm,
+ dimScaledLocalDistXY=[1.8, 1.8]) # dimension scaled dist is good for hybrid sectors with different module size
+algorithms.append(ce_ecal_cl)
+
+ce_ecal_clreco = RecoCoG("ce_ecal_clreco",
+ inputHitCollection=ce_ecal_cl.inputHitCollection,
+ inputProtoClusterCollection=ce_ecal_cl.outputProtoClusterCollection,
+ outputClusterCollection="EcalEndcapNClusters",
+ mcHits="EcalEndcapNHits",
+ samplingFraction=0.998, # this accounts for a small fraction of leakage
+ logWeightBase=4.6)
+algorithms.append(ce_ecal_clreco)
+
+# Endcap Sampling Ecal
+ci_ecal_daq = dict(
+ dynamicRangeADC=50.*units.MeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=10)
+
+ci_ecal_digi = CalHitDigi("ci_ecal_digi",
+ inputHitCollection="EcalEndcapPHits",
+ outputHitCollection="EcalEndcapPRawHits",
+ **ci_ecal_daq)
+algorithms.append(ci_ecal_digi)
+
+ci_ecal_reco = CalHitReco("ci_ecal_reco",
+ inputHitCollection=ci_ecal_digi.outputHitCollection,
+ outputHitCollection="EcalEndcapPRecHits",
+ thresholdFactor=5.0,
+ **ci_ecal_daq)
+algorithms.append(ci_ecal_reco)
+
+# merge hits in different layer (projection to local x-y plane)
+ci_ecal_merger = CalHitsMerger("ci_ecal_merger",
+ inputHitCollection=ci_ecal_reco.outputHitCollection,
+ outputHitCollection="EcalEndcapPRecMergedHits",
+ fields=["fiber_x", "fiber_y"],
+ fieldRefNumbers=[1, 1],
+ # fields=["layer", "slice"],
+ # fieldRefNumbers=[1, 0],
+ readoutClass="EcalEndcapPHits")
+algorithms.append(ci_ecal_merger)
+
+ci_ecal_cl = IslandCluster("ci_ecal_cl",
+ inputHitCollection=ci_ecal_merger.outputHitCollection,
+ outputProtoClusterCollection="EcalEndcapPProtoClusters",
+ splitCluster=False,
+ minClusterCenterEdep=10.*units.MeV,
+ localDistXY=[10*units.mm, 10*units.mm])
+algorithms.append(ci_ecal_cl)
+
+ci_ecal_clreco = RecoCoG("ci_ecal_clreco",
+ inputHitCollection=ci_ecal_cl.inputHitCollection,
+ inputProtoClusterCollection=ci_ecal_cl.outputProtoClusterCollection,
+ outputClusterCollection="EcalEndcapPClusters",
+ mcHits="EcalEndcapPHits",
+ logWeightBase=6.2,
+ samplingFraction=ci_ecal_sf)
+algorithms.append(ci_ecal_clreco)
+
+# Output
+podout = PodioOutput("out", filename=output_rec)
+podout.outputCommands = [
+ "keep *",
+ "drop *Hits",
+ "keep *RecHits",
+ "keep *Layers",
+ "keep *Clusters",
+ "drop *ProtoClusters",
+ "drop outputParticles",
+ "drop InitTrackParams",
+ ] + [ "drop " + c for c in sim_coll]
+algorithms.append(podout)
+
+ApplicationMgr(
+ TopAlg = algorithms,
+ EvtSel = 'NONE',
+ EvtMax = n_events,
+ ExtSvc = services,
+ OutputLevel = WARNING,
+ AuditAlgorithms = True
+ )
diff --git a/options/reconstruction.hcal.py b/options/reconstruction.hcal.py
new file mode 100644
index 00000000..0bf946c6
--- /dev/null
+++ b/options/reconstruction.hcal.py
@@ -0,0 +1,170 @@
+from Gaudi.Configuration import *
+
+from Configurables import ApplicationMgr, AuditorSvc, EICDataSvc, PodioOutput, GeoSvc
+from GaudiKernel import SystemOfUnits as units
+from GaudiKernel.SystemOfUnits import MeV, GeV, mm, cm, mrad
+
+import json
+
+detector_name = "athena"
+if "JUGGLER_DETECTOR" in os.environ :
+ detector_name = str(os.environ["JUGGLER_DETECTOR"])
+
+detector_path = ""
+if "DETECTOR_PATH" in os.environ :
+ detector_path = str(os.environ["DETECTOR_PATH"])
+
+compact_path = os.path.join(detector_path, detector_name)
+
+# CAL reconstruction
+# get sampling fractions from system environment variable
+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_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"])
+
+# services
+services = []
+# geometry service
+services.append(GeoSvc("GeoSvc", detectors=["{}.xml".format(compact_path)], OutputLevel=WARNING))
+# data service
+services.append(EICDataSvc("EventDataSvc", inputs=input_sims, OutputLevel=WARNING))
+
+# juggler components
+from Configurables import PodioInput
+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",
+ "HcalEndcapPHits",
+ "HcalEndcapNHits",
+]
+
+# list of algorithms
+algorithms = []
+
+# input
+podin = PodioInput("PodioReader", collections=sim_coll)
+algorithms.append(podin)
+
+# Hcal Hadron Endcap
+ci_hcal_daq = dict(
+ dynamicRangeADC=50.*units.MeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=10)
+
+ci_hcal_digi = CalHitDigi("ci_hcal_digi",
+ inputHitCollection="HcalEndcapPHits",
+ outputHitCollection="HcalEndcapPRawHits",
+ **ci_hcal_daq)
+algorithms.append(ci_hcal_digi)
+
+ci_hcal_reco = CalHitReco("ci_hcal_reco",
+ inputHitCollection=ci_hcal_digi.outputHitCollection,
+ outputHitCollection="HcalEndcapPRecHits",
+ thresholdFactor=5.0,
+ **ci_hcal_daq)
+algorithms.append(ci_hcal_reco)
+
+ci_hcal_merger = CalHitsMerger("ci_hcal_merger",
+ inputHitCollection=ci_hcal_reco.outputHitCollection,
+ outputHitCollection="HcalEndcapPMergedHits",
+ readoutClass="HcalEndcapPHits",
+ fields=["layer", "slice"],
+ fieldRefNumbers=[1, 0])
+algorithms.append(ci_hcal_merger)
+
+ci_hcal_cl = IslandCluster("ci_hcal_cl",
+ inputHitCollection=ci_hcal_merger.outputHitCollection,
+ outputProtoClusterCollection="HcalEndcapPProtoClusters",
+ splitCluster=False,
+ minClusterCenterEdep=30.*units.MeV,
+ localDistXY=[15.*units.cm, 15.*units.cm])
+algorithms.append(ci_hcal_cl)
+
+ci_hcal_clreco = RecoCoG("ci_hcal_clreco",
+ inputHitCollection=ci_hcal_cl.inputHitCollection,
+ inputProtoClusterCollection=ci_hcal_cl.outputProtoClusterCollection,
+ outputClusterCollection="HcalEndcapPClusters",
+ mcHits="HcalEndcapPHits",
+ logWeightBase=6.2,
+ samplingFraction=ci_hcal_sf)
+algorithms.append(ci_hcal_clreco)
+
+# Hcal Electron Endcap
+ce_hcal_daq = dict(
+ dynamicRangeADC=50.*units.MeV,
+ capacityADC=32768,
+ pedestalMean=400,
+ pedestalSigma=10)
+
+ce_hcal_digi = CalHitDigi("ce_hcal_digi",
+ inputHitCollection="HcalEndcapNHits",
+ outputHitCollection="HcalEndcapNRawHits",
+ **ce_hcal_daq)
+algorithms.append(ce_hcal_digi)
+
+ce_hcal_reco = CalHitReco("ce_hcal_reco",
+ inputHitCollection=ce_hcal_digi.outputHitCollection,
+ outputHitCollection="HcalEndcapNRecHits",
+ thresholdFactor=5.0,
+ **ce_hcal_daq)
+algorithms.append(ce_hcal_reco)
+
+ce_hcal_merger = CalHitsMerger("ce_hcal_merger",
+ inputHitCollection=ce_hcal_reco.outputHitCollection,
+ outputHitCollection="HcalEndcapNMergedHits",
+ readoutClass="HcalEndcapNHits",
+ fields=["layer", "slice"],
+ fieldRefNumbers=[1, 0])
+algorithms.append(ce_hcal_merger)
+
+ce_hcal_cl = IslandCluster("ce_hcal_cl",
+ inputHitCollection=ce_hcal_merger.outputHitCollection,
+ outputProtoClusterCollection="HcalEndcapNProtoClusters",
+ splitCluster=False,
+ minClusterCenterEdep=30.*units.MeV,
+ localDistXY=[15.*units.cm, 15.*units.cm])
+algorithms.append(ce_hcal_cl)
+
+ce_hcal_clreco = RecoCoG("ce_hcal_clreco",
+ inputHitCollection=ce_hcal_cl.inputHitCollection,
+ inputProtoClusterCollection=ce_hcal_cl.outputProtoClusterCollection,
+ outputClusterCollection="HcalEndcapNClusters",
+ mcHits="HcalEndcapNHits",
+ logWeightBase=6.2,
+ samplingFraction=ce_hcal_sf)
+algorithms.append(ce_hcal_clreco)
+
+# Output
+podout = PodioOutput("out", filename=output_rec)
+podout.outputCommands = [
+ "keep *",
+ "drop *Hits",
+ "keep *RecHits",
+ "keep *Layers",
+ "keep *Clusters",
+ "drop *ProtoClusters",
+ "drop outputParticles",
+ "drop InitTrackParams",
+ ] + [ "drop " + c for c in sim_coll]
+algorithms.append(podout)
+
+ApplicationMgr(
+ TopAlg = algorithms,
+ EvtSel = 'NONE',
+ EvtMax = n_events,
+ ExtSvc = services,
+ OutputLevel = WARNING,
+ AuditAlgorithms = True
+ )
--
GitLab