diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 3c94364965749e8867b0de875c30afb2e96b8397..eeb2c0e7671ec69732bf55e4e9860c3a1f02fd18 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -45,7 +45,6 @@ common:detector: - print_env.sh .rec_benchmark: - image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG needs: - ["common:detector"] before_script: diff --git a/benchmarks/clustering/full_cal_clusters.sh b/benchmarks/clustering/full_cal_clusters.sh index 71b4db4a029d51b1d973bfdaeae87132c8823dd7..ba97c36e2507e198c081919465cd71d6ace364ca 100644 --- a/benchmarks/clustering/full_cal_clusters.sh +++ b/benchmarks/clustering/full_cal_clusters.sh @@ -122,23 +122,14 @@ FULL_CAL_OPTION_DIR=benchmarks/clustering/options FULL_CAL_SCRIPT_DIR=benchmarks/clustering/scripts # Run Juggler -# Digitization +# Digitization & Clustering xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \ - gaudirun.py ${FULL_CAL_OPTION_DIR}/full_cal_digi.py + gaudirun.py ${FULL_CAL_OPTION_DIR}/full_cal_reco.py if [[ "$?" -ne "0" ]] ; then echo "ERROR running digitization (juggler)" exit 1 fi -# 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 diff --git a/benchmarks/clustering/options/calorimeter_clustering.py b/benchmarks/clustering/options/deprecated/calorimeter_clustering.py similarity index 91% rename from benchmarks/clustering/options/calorimeter_clustering.py rename to benchmarks/clustering/options/deprecated/calorimeter_clustering.py index f2ac4f2599bbf5fc57246da75bd47ce0878bb151..d4431f163fec155d2f19b0e64a760963375c6211 100644 --- a/benchmarks/clustering/options/calorimeter_clustering.py +++ b/benchmarks/clustering/options/deprecated/calorimeter_clustering.py @@ -22,7 +22,7 @@ input_sim_file = str(os.environ["JUGGLER_SIM_FILE"]) output_rec_file = str(os.environ["JUGGLER_REC_FILE"]) n_events = str(os.environ["JUGGLER_N_EVENTS"]) -geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(detector_name)]) +geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(detector_name)], OutputLevel=INFO) podioevent = EICDataSvc("EventDataSvc", inputs=[input_sim_file]) from Configurables import PodioInput @@ -70,7 +70,7 @@ ecal_reco = EMCalReconstruction("ecal_reco", ec_barrel_cluster = IslandCluster("ec_barrel_cluster", inputHitCollection="RecEcalBarrelHits", - outputClusterCollection="EcalBarrelClusters", + outputClusterCollection="EcalBarrelProtoClusters", splitHitCollection="splitEcalBarrelHitCollection", minClusterCenterEdep=1*units.MeV, groupRange=2.0, @@ -78,7 +78,8 @@ ec_barrel_cluster = IslandCluster("ec_barrel_cluster", crystal_ec_cluster = IslandCluster("crystal_ec_cluster", inputHitCollection="RecoEcalHits", - outputClusterCollection="EcalClusters", + outputClusterCollection="EcalProtoClusters", + splitHitCollection="splitEcalHitCollection", minClusterCenterEdep=30*units.MeV, groupRange=2.0, OutputLevel=DEBUG) @@ -89,12 +90,14 @@ simple_cluster = SimpleClustering("simple_cluster", OutputLevel=DEBUG) ec_barrel_clusterreco = RecoCoG("ec_barrel_clusterreco", - clusterCollection="EcalBarrelClusters", + inputClusterCollection="EcalBarrelProtoClusters", + outputClusterCollection="EcalBarrelClusters", logWeightBase=6.2, samplingFraction=sf) clusterreco = RecoCoG("cluster_reco", - clusterCollection="EcalClusters", + inputClusterCollection="EcalProtoClusters", + outputClusterCollection="EcalClusters", logWeightBase=4.2, moduleDimZName="CrystalBox_z_length", samplingFraction=sf, diff --git a/benchmarks/clustering/options/full_cal_clusters.py b/benchmarks/clustering/options/deprecated/full_cal_clusters.py similarity index 86% rename from benchmarks/clustering/options/full_cal_clusters.py rename to benchmarks/clustering/options/deprecated/full_cal_clusters.py index 0b6b34830d43c35e5b2ade2de011b0f9e83d4444..c0918e011518fbb7dd1780e357e8535aa23320cc 100644 --- a/benchmarks/clustering/options/full_cal_clusters.py +++ b/benchmarks/clustering/options/deprecated/full_cal_clusters.py @@ -27,7 +27,7 @@ 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)]) +geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(compact_path)], OutputLevel=INFO) # data service podioevent = EICDataSvc("EventDataSvc", inputs=input_digis) @@ -67,14 +67,15 @@ podout = PodioOutput("out", filename=output_rec) ce_ecal_cl = IslandCluster("ce_ecal_cl", # OutputLevel=DEBUG, inputHitCollection="EcalEndcapNHitsReco", - outputClusterCollection="EcalEndcapNClusters", + outputClusterCollection="EcalEndcapNProtoClusters", splitHitCollection="EcalEndcapNHitsSplit", splitCluster=False, minClusterCenterEdep=30*MeV, groupRanges=[2.2*cm, 2.2*cm]) ce_ecal_clreco = RecoCoG("ce_ecal_clreco", - clusterCollection="EcalEndcapNClusters", + inputClusterCollection="EcalEndcapNProtoClusters", + outputClusterCollection="EcalEndcapNClusters", samplingFraction=0.998, # this accounts for a small fraction of leakage logWeightBase=4.6) @@ -90,14 +91,15 @@ ci_ecal_merger = CalHitsMerger("ci_ecal_merger", ci_ecal_cl = IslandCluster("ci_ecal_cl", inputHitCollection="EcalEndcapPHitsRecoXY", - outputClusterCollection="EcalEndcapPClusters", + outputClusterCollection="EcalEndcapPProtoClusters", splitHitCollection="EcalEndcapPHitsSplit", splitCluster=False, minClusterCenterEdep=30.*MeV, groupRanges=[5*mm, 5*mm]) ci_ecal_clreco = RecoCoG("ci_ecal_clreco", - clusterCollection="EcalEndcapPClusters", + inputClusterCollection="EcalEndcapPProtoClusters", + outputClusterCollection="EcalEndcapPClusters", logWeightBase=6.2, samplingFraction=ce_ecal_sf) @@ -105,14 +107,16 @@ ci_ecal_clreco = RecoCoG("ci_ecal_clreco", # Central Barrel Ecal (Imaging Cal.) cb_ecal_cl = ImagingCluster("cb_ecal_cl", inputHitCollection="EcalBarrelHitsReco", - outputClusterCollection="EcalBarrelClusters", + outputClusterCollection="EcalBarrelProtoClusters", + 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 cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco", samplingFraction=cb_ecal_sf, - inputClusterCollection="EcalBarrelClusters", + inputClusterCollection="EcalBarrelProtoClusters", + outputClusterCollection="EcalBarrelClusters", outputLayerCollection="EcalBarrelLayers") @@ -127,14 +131,15 @@ cb_hcal_merger = CalHitsMerger("cb_hcal_merger", cb_hcal_cl = IslandCluster("cb_hcal_cl", inputHitCollection="HcalBarrelHitsRecoXY", - outputClusterCollection="HcalBarrelClusters", + outputClusterCollection="HcalBarrelProtoClusters", splitHitCollection="HcalBarrelHitsSplit", splitCluster=False, minClusterCenterEdep=30.*MeV, groupRanges=[15.*cm, 15.*cm]) cb_hcal_clreco = RecoCoG("cb_hcal_clreco", - clusterCollection="HcalBarrelClusters", + inputClusterCollection="HcalBarrelProtoClusters", + outputClusterCollection="HcalBarrelClusters", logWeightBase=6.2, samplingFraction=cb_hcal_sf) @@ -149,14 +154,15 @@ ci_hcal_merger = CalHitsMerger("ci_hcal_merger", ci_hcal_cl = IslandCluster("ci_hcal_cl", inputHitCollection="HcalHadronEndcapHitsRecoXY", - outputClusterCollection="HcalHadronEndcapClusters", + outputClusterCollection="HcalHadronEndcapProtoClusters", splitHitCollection="HcalHadronEndcapHitsSplit", splitCluster=False, minClusterCenterEdep=30.*MeV, groupRanges=[15.*cm, 15.*cm]) ci_hcal_clreco = RecoCoG("ci_hcal_clreco", - clusterCollection="HcalHadronEndcapClusters", + inputClusterCollection="HcalHadronEndcapProtoClusters", + outputClusterCollection="HcalHadronEndcapClusters", logWeightBase=6.2, samplingFraction=ci_hcal_sf) @@ -170,14 +176,15 @@ ce_hcal_merger = CalHitsMerger("ce_hcal_merger", ce_hcal_cl = IslandCluster("ce_hcal_cl", inputHitCollection="HcalElectronEndcapHitsRecoXY", - outputClusterCollection="HcalElectronEndcapClusters", + outputClusterCollection="HcalElectronEndcapProtoClusters", splitHitCollection="HcalElectronEndcapHitsSplit", splitCluster=False, minClusterCenterEdep=30.*MeV, groupRanges=[15.*cm, 15.*cm]) ce_hcal_clreco = RecoCoG("ce_hcal_clreco", - clusterCollection="HcalElectronEndcapClusters", + inputClusterCollection="HcalElectronEndcapProtoClusters", + outputClusterCollection="HcalElectronEndcapClusters", logWeightBase=6.2, samplingFraction=ce_hcal_sf) diff --git a/benchmarks/clustering/options/full_cal_digi.py b/benchmarks/clustering/options/deprecated/full_cal_digi.py similarity index 99% rename from benchmarks/clustering/options/full_cal_digi.py rename to benchmarks/clustering/options/deprecated/full_cal_digi.py index 8850e17c7fbdd2513a5adf3555dc14f1d5281563..347e14e246dd8d7e869188f388c551b28defe07a 100644 --- a/benchmarks/clustering/options/full_cal_digi.py +++ b/benchmarks/clustering/options/deprecated/full_cal_digi.py @@ -19,7 +19,7 @@ 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)]) +geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(compact_path)], OutputLevel=INFO) # data service podioevent = EICDataSvc("EventDataSvc", inputs=input_sims) diff --git a/benchmarks/clustering/options/hcal_clustering.py b/benchmarks/clustering/options/deprecated/hcal_clustering.py similarity index 100% rename from benchmarks/clustering/options/hcal_clustering.py rename to benchmarks/clustering/options/deprecated/hcal_clustering.py diff --git a/benchmarks/clustering/options/full_cal_reco.py b/benchmarks/clustering/options/full_cal_reco.py new file mode 100644 index 0000000000000000000000000000000000000000..1c7d0508e61963abbd5f535a95c8e942c2b390e3 --- /dev/null +++ b/benchmarks/clustering/options/full_cal_reco.py @@ -0,0 +1,309 @@ +''' + 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("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)) +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_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=["{}.xml".format(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__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", + "EcalEndcapNHits", + "EcalEndcapPHits", + "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="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_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 + readoutClass="EcalEndcapNHits", + sectorField="sector") + +ce_ecal_cl = IslandCluster("ce_ecal_cl", + # OutputLevel=DEBUG, + inputHitCollection="EcalEndcapNHitsReco", + outputHitCollection="EcalEndcapNClusterHits", + splitCluster=False, + minClusterCenterEdep=30*MeV, + groupRanges=[2.2*cm, 2.2*cm]) + +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) + +# Endcap Sampling Ecal +ci_ecal_digi = CalHitDigi("ci_ecal_digi", + inputHitCollection="EcalEndcapPHits", + outputHitCollection="EcalEndcapPHitsDigi", + dynamicRangeADC=50.*MeV, + capacityADC=32768, + pedestalMean=400, + pedestalSigma=10) + +ci_ecal_reco = CalHitReco("ci_ecal_reco", + inputHitCollection="EcalEndcapPHitsDigi", + outputHitCollection="EcalEndcapPHitsReco", + dynamicRangeADC=50.*MeV, + capacityADC=32768, + pedestalMean=400, + pedestalSigma=10, + thresholdFactor=5.0) +# merge hits in different layer (projection to local x-y plane) +ci_ecal_merger = CalHitsMerger("ci_ecal_merger", + inputHitCollection="EcalEndcapPHitsReco", + outputHitCollection="EcalEndcapPHitsRecoXY", + fields=["layer", "slice"], + fieldRefNumbers=[1, 0], + readoutClass="EcalEndcapPHits") + +ci_ecal_cl = IslandCluster("ci_ecal_cl", + inputHitCollection="EcalEndcapPHitsRecoXY", + outputHitCollection="EcalEndcapPClusterHits", + splitCluster=False, + minClusterCenterEdep=30.*MeV, + groupRanges=[5*mm, 5*mm]) + +ci_ecal_clreco = RecoCoG("ci_ecal_clreco", + inputHitCollection="EcalEndcapPClusterHits", + outputClusterCollection="EcalEndcapPClusters", + 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", + 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 +cb_ecal_clreco = ImagingClusterReco("cb_ecal_clreco", + samplingFraction=cb_ecal_sf, + inputHitCollection="EcalBarrelClusterHits", + outputClusterCollection="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, + readoutClass="HcalBarrelHits", + layerField="layer", + sectorField="module") +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", + outputHitCollection="HcalBarrelClusterHits", + splitCluster=False, + minClusterCenterEdep=30.*MeV, + groupRanges=[15.*cm, 15.*cm]) + +cb_hcal_clreco = RecoCoG("cb_hcal_clreco", + inputHitCollection="HcalBarrelClusterHits", + outputClusterCollection="HcalBarrelClusters", + logWeightBase=6.2, + samplingFraction=cb_hcal_sf) + + +# 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) +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", + outputHitCollection="HcalHadronEndcapClusterHits", + splitCluster=False, + minClusterCenterEdep=30.*MeV, + groupRanges=[15.*cm, 15.*cm]) + +ci_hcal_clreco = RecoCoG("ci_hcal_clreco", + inputHitCollection="HcalHadronEndcapClusterHits", + outputClusterCollection="HcalHadronEndcapClusters", + logWeightBase=6.2, + samplingFraction=ci_hcal_sf) + +# 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) +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", + outputHitCollection="HcalElectronEndcapClusterHits", + splitCluster=False, + minClusterCenterEdep=30.*MeV, + groupRanges=[15.*cm, 15.*cm]) + +ce_hcal_clreco = RecoCoG("ce_hcal_clreco", + inputHitCollection="HcalElectronEndcapClusterHits", + outputClusterCollection="HcalElectronEndcapClusters", + logWeightBase=6.2, + samplingFraction=ce_hcal_sf) + + +podout.outputCommands = ['drop *', 'keep mcparticles*', 'keep *Reco', 'keep *Digi', 'keep *Clusters', 'keep *Layers', 'drop *ProtoClusters'] + +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, + podout], + EvtSel = 'NONE', + EvtMax = n_events, + ExtSvc = [podioevent], + OutputLevel=DEBUG +) diff --git a/benchmarks/ecal/options/crystal_calorimeter_reco.py b/benchmarks/ecal/options/crystal_calorimeter_reco.py index e9adbe716450acb5187930921f9262c9d8362506..30b563d5ffa7d341dc5eacb6be6f3692f366d724 100644 --- a/benchmarks/ecal/options/crystal_calorimeter_reco.py +++ b/benchmarks/ecal/options/crystal_calorimeter_reco.py @@ -31,7 +31,7 @@ n_events = 100 if "JUGGLER_N_EVENTS" in os.environ : n_events = str(os.environ["JUGGLER_N_EVENTS"]) -geo_service = GeoSvc("GeoSvc", detectors=["{}/{}.xml".format(detector_path,detector_name)]) +geo_service = GeoSvc("GeoSvc", detectors=["{}/{}.xml".format(detector_path,detector_name)], OutputLevel=INFO) podioevent = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=INFO) from Configurables import PodioInput @@ -62,12 +62,13 @@ emcalreco = CrystalEndcapsReco("ecal_reco", emcalcluster = IslandCluster("emcal_cluster", inputHitCollection="RecoEcalHits", - outputClusterCollection="EcalClusters", + outputHitCollection="EcalClusterHits", minClusterCenterEdep=50.0*units.MeV, groupRange=2.0) clusterreco = RecoCoG("cluster_reco", - clusterCollection="EcalClusters", + inputHitCollection="EcalClusterHits", + outputClusterCollection="EcalClusters", logWeightBase=4.2, moduleDimZName="CrystalModule_sz") diff --git a/benchmarks/ecal/options/example_crystal.py b/benchmarks/ecal/options/deprecated/example_crystal.py similarity index 81% rename from benchmarks/ecal/options/example_crystal.py rename to benchmarks/ecal/options/deprecated/example_crystal.py index 5c41b001183458d95a582a159d301ef0f25d137b..7f9e6ee4c008e7ff942b17b0c101067b9c963eb1 100644 --- a/benchmarks/ecal/options/example_crystal.py +++ b/benchmarks/ecal/options/deprecated/example_crystal.py @@ -25,11 +25,9 @@ from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG podioinput = PodioInput("PodioReader", collections=["mcparticles","CrystalEcalHits"], OutputLevel=DEBUG) emcaldigi = CrystalEndcapsDigi("ecal_digi", inputHitCollection="CrystalEcalHits", outputHitCollection="RawDigiEcalHits") -emcalreco = CrystalEndcapsReco("ecal_reco", inputHitCollection="RawDigiEcalHits", outputHitCollection="RecoEcalHits", - minModuleEdep=1.0*units.MeV) -emcalcluster = IslandCluster("emcal_cluster", inputHitCollection="RecoEcalHits", outputClusterCollection="EcalClusters", - minClusterCenterEdep=30*units.MeV, groupRange=2.0) -clusterreco = RecoCoG("cluster_reco", clusterCollection="EcalClusters", logWeightBase=4.2, moduleDimZName="CrystalModule_sz") +emcalreco = CrystalEndcapsReco("ecal_reco", inputHitCollection="RawDigiEcalHits", outputHitCollection="RecoEcalHits", minModuleEdep=1.0*units.MeV) +emcalcluster = IslandCluster("emcal_cluster", inputHitCollection="RecoEcalHits", outputHitCollection="EcalClusterHits", minClusterCenterEdep=30*units.MeV, groupRange=2.0) +clusterreco = RecoCoG("cluster_reco", inputGitCollection="EcalClusterHits", outputClusterCollection="EcalClusters", logWeightBase=4.2, moduleDimZName="CrystalModule_sz") out = PodioOutput("out", filename=output_rec_file) diff --git a/benchmarks/ecal/options/emcal_barrel_reco.py b/benchmarks/ecal/options/emcal_barrel_reco.py index 2abf07695ba4e30a87395af0e70cdafd4c296375..964e5c3797b2a0b670b6ca4cd4c57709b9147554 100644 --- a/benchmarks/ecal/options/emcal_barrel_reco.py +++ b/benchmarks/ecal/options/emcal_barrel_reco.py @@ -38,7 +38,7 @@ sf = 1.0 if "CB_EMCAL_SAMP_FRAC" in os.environ : sf = str(os.environ["CB_EMCAL_SAMP_FRAC"]) -geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(detector_name)]) +geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(detector_name)], OutputLevel=INFO) podioevent = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=DEBUG) from Configurables import PodioInput @@ -99,15 +99,18 @@ embarrelzmerger = CalorimeterHitsMerger("ecal_barrel_z_merger", # Clustering embarrelcluster = IslandCluster("ecal_barrel_cluster", inputHitCollection="RecoEcalBarrelHitsXY", - outputClusterCollection="EcalBarrelClusters", + outputHitCollection="EcalBarrelClusterHits", minClusterCenterEdep=0.5*units.MeV, splitCluster=False, - groupRanges=[2.0*units.cm, 2.0*units.cm, 2.0*units.cm]) + groupRanges=[2.0*units.cm, 2.0*units.cm, 2.0*units.cm], + OutputLevel=DEBUG) # Reconstruct the cluster with Center of Gravity method embarrelclusterreco = RecoCoG("ecal_barrel_clusterreco", - clusterCollection="EcalBarrelClusters", + inputHitCollection="EcalBarrelClusterHits", + outputClusterCollection="EcalBarrelClusters", logWeightBase=6.2, - samplingFraction=sf) + samplingFraction=sf, + OutputLevel=DEBUG) out = PodioOutput("out", filename=output_rec_file) diff --git a/benchmarks/ecal/options/full_em_calorimeter_reco.py b/benchmarks/ecal/options/full_em_calorimeter_reco.py index 8bbf1bdfcb13151d013ac6772cd4a24cbae68d45..67c531fd7940606e06fdfb571f10800c8ac8fcda 100644 --- a/benchmarks/ecal/options/full_em_calorimeter_reco.py +++ b/benchmarks/ecal/options/full_em_calorimeter_reco.py @@ -28,7 +28,7 @@ n_events = 100 if "JUGGLER_N_EVENTS" in os.environ : n_events = str(os.environ["JUGGLER_N_EVENTS"]) -geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(detector_name)]) +geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(detector_name)], OutputLevel=INFO) podioevent = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=DEBUG) from Configurables import PodioInput @@ -87,12 +87,13 @@ emcalreco = CrystalEndcapsReco("ecal_reco", emcalcluster = IslandCluster("emcal_cluster", inputHitCollection="RecoEcalHits", - outputClusterCollection="EcalClusters", + outputHitCollection="EcalClusterHits", minClusterCenterEdep=50.0*units.MeV, groupRange=2.0) clusterreco = RecoCoG("cluster_reco", - clusterCollection="EcalClusters", + inputHitCollection="EcalClusterHits", + outputClusterCollection="EcalClusters", logWeightBase=4.2, moduleDimZName="CrystalModule_sz") diff --git a/benchmarks/ecal/scripts/emcal_barrel_electrons_analysis.cxx b/benchmarks/ecal/scripts/emcal_barrel_electrons_analysis.cxx index 7edea3163bb4a8d2abeec49ba62267b27398380b..bab641eb0b9bb6b4de699446962d54751a27f5d2 100644 --- a/benchmarks/ecal/scripts/emcal_barrel_electrons_analysis.cxx +++ b/benchmarks/ecal/scripts/emcal_barrel_electrons_analysis.cxx @@ -6,21 +6,21 @@ #include "ROOT/RDataFrame.hxx" #include <iostream> -#include "dd4pod/Geant4ParticleCollection.h" #include "dd4pod/CalorimeterHitCollection.h" -#include "eicd/RawCalorimeterHitCollection.h" -#include "eicd/RawCalorimeterHitData.h" +#include "dd4pod/Geant4ParticleCollection.h" #include "eicd/CalorimeterHitCollection.h" #include "eicd/CalorimeterHitData.h" #include "eicd/ClusterCollection.h" #include "eicd/ClusterData.h" +#include "eicd/RawCalorimeterHitCollection.h" +#include "eicd/RawCalorimeterHitData.h" #include "TCanvas.h" -#include "TStyle.h" -#include "TMath.h" -#include "TH1.h" #include "TF1.h" +#include "TH1.h" #include "TH1D.h" +#include "TMath.h" +#include "TStyle.h" using ROOT::RDataFrame; using namespace ROOT::VecOps; @@ -44,37 +44,38 @@ void emcal_barrel_electrons_analysis(const char* input_fname = "rec_emcal_barrel // Thrown Energy [GeV] auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) { std::vector<double> result; - result.push_back(TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass)); - return result; + result.push_back(TMath::Sqrt(input[2].psx * input[2].psx + input[2].psy * input[2].psy + + input[2].psz * input[2].psz + input[2].mass * input[2].mass)); + return result; }; // Reconstructed Energy [GeV] in XY merger - auto ErecXY = [] (const std::vector<eic::CalorimeterHitData> & evt) { + auto ErecXY = [](const std::vector<eic::CalorimeterHitData>& evt) { std::vector<double> result; - auto total_eng = 0.0; - for (const auto& i: evt) + auto total_eng = 0.0; + for (const auto& i : evt) total_eng += i.energy; result.push_back(total_eng / 1.e+3); return result; }; // Reconstructed Energy [GeV] in Z merger - auto ErecZ = [] (const std::vector<eic::CalorimeterHitData> & evt) { + auto ErecZ = [](const std::vector<eic::CalorimeterHitData>& evt) { std::vector<double> result; - auto total_eng = 0.0; - for (const auto& i: evt) + auto total_eng = 0.0; + for (const auto& i : evt) total_eng += i.energy; result.push_back(total_eng / 1.e+3); return result; }; // Number of Clusters - auto ncluster = [] (const std::vector<eic::ClusterData>& evt) {return (int) evt.size(); }; + auto ncluster = [](const std::vector<eic::ClusterData>& evt) { return (int)evt.size(); }; // Cluster Energy [GeV] - auto Ecluster = [] (const std::vector<eic::ClusterData>& evt) { + auto Ecluster = [](const std::vector<eic::ClusterData>& evt) { std::vector<double> result; - for (const auto& i: evt) + for (const auto& i : evt) result.push_back(i.energy / 1.e+3); return result; }; @@ -84,34 +85,35 @@ void emcal_barrel_electrons_analysis(const char* input_fname = "rec_emcal_barrel std::vector<double> result; for (const auto& E1 : thrown) { for (const auto& E2 : sampled) - result.push_back(E2/E1); + result.push_back(E2 / E1); } return result; }; // Define variables - auto d1 = d0.Define("Ethr", Ethr, {"mcparticles2"}) - .Define("ErecXY", ErecXY, {"RecoEcalBarrelHitsXY"}) - .Define("ErecZ", ErecZ, {"RecoEcalBarrelHitsZ"}) - .Define("ncluster", ncluster, {"EcalBarrelClusters"}) - .Define("Ecluster", Ecluster, {"EcalBarrelClusters"}) - .Define("fsam", fsam, {"Ecluster","Ethr"}) - ; + auto d1 = d0.Define("Ethr", Ethr, {"mcparticles2"}) + .Define("ErecXY", ErecXY, {"RecoEcalBarrelHitsXY"}) + .Define("ErecZ", ErecZ, {"RecoEcalBarrelHitsZ"}) + .Define("ncluster", ncluster, {"EcalBarrelClusters"}) + .Define("Ecluster", Ecluster, {"EcalBarrelClusters"}) + .Define("fsam", fsam, {"Ecluster", "Ethr"}); // Define Histograms - auto hEthr = d1.Histo1D({"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 6.5}, "Ethr"); - auto hErecXY = d1.Histo1D({"hErecXY", "Reconstructed Energy in XY merger; Reconstructed Energy [GeV]; Events", 100, 0.0, 6.5}, "ErecXY"); - auto hErecZ = d1.Histo1D({"hErecZ", "Reconstructed Energy in Z merger; Reconstructed Energy [GeV]; Events", 100, 0.0, 6.5}, "ErecZ"); - auto hNCluster = d1.Histo1D({"hNCluster", "Number of Clusters; # of Clusters; Events", 100, 0.0, 100.0}, "ncluster"); - auto hEcluster = d1.Histo1D({"hEcluster", "Cluster Energy; Cluster Energy [GeV]; Events", 100, 0.0, 6.5}, "Ecluster"); - auto hfsam = d1.Histo1D({"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1}, "fsam"); + auto hEthr = d1.Histo1D({"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 6.5}, "Ethr"); + auto hErecXY = d1.Histo1D( + {"hErecXY", "Reconstructed Energy in XY merger; Reconstructed Energy [GeV]; Events", 100, 0.0, 6.5}, "ErecXY"); + auto hErecZ = d1.Histo1D( + {"hErecZ", "Reconstructed Energy in Z merger; Reconstructed Energy [GeV]; Events", 100, 0.0, 6.5}, "ErecZ"); + auto hNCluster = d1.Histo1D({"hNCluster", "Number of Clusters; # of Clusters; Events", 100, 0.0, 100.0}, "ncluster"); + auto hEcluster = d1.Histo1D({"hEcluster", "Cluster Energy; Cluster Energy [GeV]; Events", 100, 0.0, 6.5}, "Ecluster"); + auto hfsam = d1.Histo1D({"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1}, "fsam"); // Event Counts - auto nevents_thrown = d1.Count(); + auto nevents_thrown = d1.Count(); std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n"; // Draw Histograms - TCanvas *c1 = new TCanvas("c1", "c1", 700, 500); + TCanvas* c1 = new TCanvas("c1", "c1", 700, 500); c1->SetLogy(1); hEthr->GetYaxis()->SetTitleOffset(1.4); hEthr->SetLineWidth(2); @@ -120,7 +122,7 @@ void emcal_barrel_electrons_analysis(const char* input_fname = "rec_emcal_barrel c1->SaveAs("results/emcal_barrel_electrons_Ethr.png"); c1->SaveAs("results/emcal_barrel_electrons_Ethr.pdf"); - TCanvas *c2 = new TCanvas("c2", "c2", 700, 500); + TCanvas* c2 = new TCanvas("c2", "c2", 700, 500); c2->SetLogy(1); hErecXY->GetYaxis()->SetTitleOffset(1.4); hErecXY->SetLineWidth(2); @@ -129,17 +131,17 @@ void emcal_barrel_electrons_analysis(const char* input_fname = "rec_emcal_barrel c2->SaveAs("results/emcal_barrel_electrons_ErecXY.png"); c2->SaveAs("results/emcal_barrel_electrons_ErecXY.pdf"); - TCanvas *c3 = new TCanvas("c3", "c3", 700, 500); + TCanvas* c3 = new TCanvas("c3", "c3", 700, 500); c3->SetLogy(1); hErecZ->GetYaxis()->SetTitleOffset(1.4); hErecZ->SetLineWidth(2); hErecZ->SetLineColor(kBlue); hErecZ->DrawClone(); - c3->SaveAs("results/emcal_electrons_ErecZ.png"); + c3->SaveAs("results/emcal_electrons_ErecZ.png"); c3->SaveAs("results/emcal_electrons_ErecZ.pdf"); - TCanvas *c4 = new TCanvas("c4", "c4", 700, 500); - c4->SetLogy(1); + TCanvas* c4 = new TCanvas("c4", "c4", 700, 500); + c4->SetLogy(1); hNCluster->GetYaxis()->SetTitleOffset(1.6); hNCluster->SetLineWidth(2); hNCluster->SetLineColor(kBlue); @@ -147,7 +149,7 @@ void emcal_barrel_electrons_analysis(const char* input_fname = "rec_emcal_barrel c4->SaveAs("results/emcal_barrel_electrons_ncluster.png"); c4->SaveAs("results/emcal_barrel_electrons_ncluster.pdf"); - TCanvas *c5 = new TCanvas("c5", "c5", 700, 500); + TCanvas* c5 = new TCanvas("c5", "c5", 700, 500); c5->SetLogy(1); hEcluster->GetYaxis()->SetTitleOffset(1.4); hEcluster->SetLineWidth(2); @@ -156,13 +158,13 @@ void emcal_barrel_electrons_analysis(const char* input_fname = "rec_emcal_barrel c5->SaveAs("results/emcal_barrel_electrons_Ecluster.png"); c5->SaveAs("results/emcal_barrel_electrons_Ecluster.pdf"); - TCanvas *c6 = new TCanvas("c6", "c6", 700, 500); + TCanvas* c6 = new TCanvas("c6", "c6", 700, 500); c6->SetLogy(1); hfsam->GetYaxis()->SetTitleOffset(1.4); hfsam->SetLineWidth(2); hfsam->SetLineColor(kBlue); - hfsam->Fit("gaus","","",0.01,0.1); - // From emcal_barrel_pions_analysis.cxx + hfsam->Fit("gaus", "", "", 0.01, 0.1); + // From emcal_barrel_pions_analysis.cxx // Commented out for now as giving issues with new container (S. Joosten) // hfsam->GetFunction("gaus")->SetLineWidth(2); // hfsam->GetFunction("gaus")->SetLineColor(kRed); diff --git a/benchmarks/imaging_ecal/options/imaging_topocluster.py b/benchmarks/imaging_ecal/options/imaging_topocluster.py index cdcbfc0238372c79c450035cfda1564532f4d880..722985117bcb120eaca5921baf673751d766cb6b 100644 --- a/benchmarks/imaging_ecal/options/imaging_topocluster.py +++ b/benchmarks/imaging_ecal/options/imaging_topocluster.py @@ -31,7 +31,7 @@ print(kwargs) # get sampling fraction from system environment variable, 1.0 by default sf = float(os.environ.get('CB_EMCAL_SAMP_FRAC', '1.0')) -geo_service = GeoSvc("GeoSvc", detectors=kwargs['compact'].split(',')) +geo_service = GeoSvc("GeoSvc", detectors=kwargs['compact'].split(','), OutputLevel=INFO) podioevent = EICDataSvc("EventDataSvc", inputs=kwargs['input'].split(','), OutputLevel=DEBUG) out = PodioOutput("out", filename=kwargs['output']) @@ -65,24 +65,24 @@ imcalreco = ImagingPixelReco("imcal_reco", sectorField="module") imcalcluster = ImagingTopoCluster("imcal_cluster", inputHitCollection="RecoEcalBarrelHits", - outputClusterCollection="EcalBarrelClusters", + outputHitCollection="EcalBarrelClusterHits", localRanges=[2.*units.mm, 2*units.mm], adjLayerRanges=[10*units.mrad, 10*units.mrad], adjLayerDiff=2, adjSectorDist=3.*units.cm) clusterreco = ImagingClusterReco("imcal_clreco", - inputClusterCollection="EcalBarrelClusters", + inputHitCollection="EcalBarrelClusterHits", outputLayerCollection="EcalBarrelClustersLayers", + outputClusterCollection="EcalBarrelClusters", samplingFraction=sf, OutputLevel=DEBUG) out.outputCommands = ["keep *"] ApplicationMgr( - TopAlg=[podioinput, copier, calcopier, imcaldigi, imcalreco, imcalcluster, out], + TopAlg=[podioinput, copier, calcopier, imcaldigi, imcalreco, imcalcluster, clusterreco, out], EvtSel='NONE', EvtMax=kwargs['nev'], ExtSvc=[podioevent], OutputLevel=DEBUG ) -